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Abstract The microscopic definition for the Cauchy stress tensor has been examined 
in the past from many different perspectives. This has led to different expressions for 
the stress tensor and consequently the "correct" definition has been a subject of debate 
and controversy. In this work, a unified framework is set up in which all existing def- 
initions can be derived, thus establishing the connections between them. The frame- 
work is based on the non-equilibrium statistical mechanics procedure introduced by 
Irving, Kirkwood and Noll, followed by spatial averaging. The Irving-Kirkwood- 
Noll procedure is extended to multi-body potentials with continuously differentiable 
extensions and generalized to non-straight bonds, which may be important for parti- 
cles with internal structure. Connections between this approach and the direct spatial 
averaging approach of Murdoch and Hardy are discussed and the Murdoch-Hardy 
procedure is systematized. Possible sources of non-uniqueness of the stress tensor, 
resulting separately from both procedures, are identified and addressed. Numerical 
experiments using molecular dynamics and lattice statics are conducted to examine 
the behavior of the resulting stress definitions including their convergence with the 
spatial averaging domain size and their symmetry properties. 
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1 Introduction 

Continuum mechanics provides an efficient theoretical framework for modeling ma- 
terials science phenomena. To characterize the behavior of materials, constitutive re- 
lations serve as an input to the continuum theory. These constitutive models have 
functional forms which must be consistent with material frame-indifference and the 
laws of thermodynamics and include parameters that are fitted to reproduce exper- 
imental observations. With the advent of modern computing power, atomistic sim- 
ulations through "numerical experiments" offer the potential for studying different 
materials and arriving at their constitutive laws from first principles. This could make 
it possible to design new materials and to improve the properties of existing mate- 
rials in a systematic fashion. To use the data obtained from an atomistic simulation 
to build a constitutive law, which is framed in the language of continuum mechan- 
ics, it is necessary to understand the connection between continuum fields and the 
underlying microscopic dynamics. 

Another arena where the connection between continuum and atomistic concepts 
is important is the field of multiscale modeling. This discipline involves the develop- 
ment of computational tools for studying problems where two or more length and/or 
time scales play a major role in determining macroscopic behavior. A prototypical 
example is fracture mechanics where the behavior of a crack is controlled by atomic- 
scale phenomena at the crack-tip, while at the same time long-range elastic stress 
fields are set up in the body. Many advances have been made in the area of multiscale 
modeling in recent years. Some common atomistic/continuum coupling methods are 
quasicontinuum 15311591 . coupling of lengthscales lf5TI . cluster quasicontinuum BP . 
bridging domain Il66l . coupled atomistics and discrete dislocations Il54l . and hetero- 
geneous multiscale methods |[T5l . to name just a few. Refer to l58l for a comparison 
of some prominent atomistic/continuum coupling multiscale methods. In a multiscale 
method, a key issue involves the transfer of information between the discrete model 
and the continuum model. It is therefore of practical interest to understand how to 
construct definitions of continuum fields for an atomistic system, to ensure a smooth 
transfer of information between the discrete and continuum domains. 

In this paper, we focus on just one aspect of the continuum-atomistic connection, 
namely the interpretation of the (Cauchy) stress tensor in a discrete system. This 
question has been explored from many different perspectives for nearly two hundred 
years and this has led to various definitions that do not appear to be consistent with 
each other. As a result, the "correct" definition for the stress tensor has been a subject 
of great debate and controversy. We begin with a brief historical survey. 

A brief history of microscopic definitions for the stress tensor 

Interest in microscopic definitions for the stress tensor dates back at least to Cauchy 
in the 1820s 06) with his aim to define stress in a crystalline solid. Cauchy's original 
definition emerges from the intuitive idea of identifying stress with the force per unit 
area carried by the bonds that cross a given surface. A comprehensive derivation 
of Cauchy's approach is given in Note B of Love's classic book on the theory of 
elasticity ll35l . Since this approach is tied to the particular surface being considered, 
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it actually constitutes a definition for the traction (or stress vector) and not for the 
stress tensor. The first definition of stress as a tensorial quantity follows from the 
works of Clausius (8) and Maxwell 140114 II in the form of virial theorem. Clausius 
states the virial theorem as 

The mean vis viva of a system is equal to its virial. 

By "vis viva" (literally "living force"), Clausius means kinetic energy, while the term 
"virial" comes from the Latin "vis" (pi. "vires") meaning force. The virial theorem 
leads to a definition for pressure in a gas. Maxwell 04OII411 extended Clausius' work 
and showed the existence of a tensorial version of the virial theorem (see Appendix 
lAl . The virial stress resulting from the virial theorem is widely used even today in 
many atomistic simulations due to its simple form and ease of computation. Unlike 
Cauchy's original definition for stress, the virial stress includes a contribution due 
to the kinetic energy of the particles. This discrepancy was addressed by Tsai l63l . 
who extended the definition given by Cauchy to finite temperature by taking into 
consideration the momentum flux passing through the surface. Let us refer to this 
stress vector as the Tsai traction. 

An alternative approach for defining the stress tensor was pioneered in the land- 
mark paper of Irving and Kirkwood ll27l . Irving and Kirkwood derived the equa- 
tions of hydrodynamics from the principles of non-equilibrium classical statistical 
mechanics and in the process established a pointwise definition for various contin- 
uum fields including the stress tensor. Although their work was indeed noteworthy, 
the stress tensor obtained involved a series expansion of the Dirac delta distribution 
which is not mathematically rigorous. Continuing their work, Noll l48l proved two 
lemmas, which allowed him to avoid the use of the Dirac delta distribution, and thus 
arrive at a closed-form expression for the stress tensor which does not involve a se- 
ries expansion. We refer to the procedure introduced by Irving and Kirkwood and 
extended by Noll as the Irving-Kirkwood-Noll procedure. Schofield and Hender- 
son Il52l highlighted the non-uniqueness present in the stress tensor derived by Irving 
and Kirkwood, and pointed out that it could result in a non-symmetric stress tensor. 
There have been several attempts to improve on the Irving and Kirkwood procedure. 
In particular, Lutsko |36| reformulated this procedure in Fourier space. An error in 
Lutsko's derivation was corrected by Cormier et al. J5)- 

Due to the stochastic nature of the Irving and Kirkwood stress, many difficulties 
arise when one tries to use their expression in atomistic simulations. To avoid these 
difficulties, Hardy and co-workers Il23lf24l and independently Murdoch 143 47l de- 
veloped a new approach that bypasses the mathematical complexity of the Irving and 
Kirkwood procedure. This is done by defining continuum fields as direct spatial av- 
erages of the discrete equations of motion using weighting functions with compact 
support. In particular, this approach leads to the so-called Hardy stress 1231 often 
used in molecular dynamics simulations. Murdoch in [45) provides an excellent de- 
scription of the spatial averaging approaches currently being used and discusses the 
non-uniqueness of the stress tensor resulting from the spatial averaging procedure. 
We refer to the direct spatial averaging approach as the Murdoch-Hardy procedure. 

Another approach, which leads to a stress tensor very similar to that obtained 
by Irving and Kirkwood is the reformulation of elasticity theory using peridynam- 
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ics ll55l . Lehoucq and Silling ll33ll have recently shown that Noll's solution is a mini- 
mum solution in a variational sense. Morante et al. |42| proposed a new approach for 
defining the stress tensor using the invariance of partition function under infinitesimal 
canonical point transformations. However, their approach is limited to equilibrium 
statistical mechanics and involves taking derivatives of delta distributions. 

We can summarize the "state of the art" for the microscopic definition of the 
stress tensor as follows. There are currently at least three definitions for the stress 
tensor which are commonly used in atomistic simulations: the virial stress, the Tsai 
traction, and the Hardy stress [68]. The importance of the Irving and Kirkwood for- 
mulation is recognized, however, it is not normally used in practice and its connection 
with the other stress definitions is not commonly understood. The difference between 
pointwise stress measures and temporal and/or spatially-averaged quantities is often 
not fully appreciated. The result is that the connection between the Cauchy stress ten- 
sor defined in continuum mechanics and its analogue, defined for a discrete system, 
remains controversial and continues to be a highly-debated problem. 

A unified framework for the microscopic definition for the stress tensor 

In this paper, a unified framework based on the Irving-Kirkwood-Noll procedure is 
established which leads to all of the major stress definitions discussed above and iden- 
tifies additional possible definitions. Since all of the definitions are obtained from a 
common framework the connections between them can be explored and analyzed and 
the uniqueness of the stress tensor can be established. An overview of the approach 
and the organization of the paper are described below. 

Before turning to the general framework, we begin in Section|2]with a derivation 
of the virial stress tensor within the framework of equilibrium statistical mechanics 
using the technique of canonical transformations. Although this derivation is quite 
different from the Irving-Kirkwood-Noll procedure, it provides insight into how the 
geometric ideas of mechanics can be used to derive the stress tensor. It also pro- 
vides a limit to which the general non-equilibrium stress tensor must converge under 
equilibrium conditions in the thermodynamic limit. This is used later to establish the 
uniqueness of the stress tensor obtained from our general unified framework. 

Next, we turn to the construction of the new unified framework. In Section [3] 
we extend the Irving-Kirkwood-Noll procedure f27] [48|, originally derived for pair 
potential interactions, to multi-body potentials. Due to the invariance of the poten- 
tial energy function with respect to the Euclidean group, it can be shown that any 
multi-body potential can be expressed as a function of distances between particles. 
When expressed in this form, we note that for a system of more than 4 particles, this 
function is only defined on a manifold since the N(N — l)/2 distances between N 
particles in M 3 are not independent for N > 5. To apply the Irving-Kirkwood-Noll 
procedure to multi-body potentials, we recognize that the potential energy function 
must be extended from its manifold to a higher-dimensional Euclidean space as a 
continuously differentiable function. We show that if such an extension exists, then 
an infinite number of equivalent extensions can be constructed using Cayley-Menger 
determinants, which describe the constraints that the distances between particles em- 
bedded in M 3 must satisfy. Then for multi-body potentials that possess continuously 
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differentiable extensions (which is the case for most practical interatomic potentials), 
we establish the key result that due to the balance of linear and angular momentum, 
the force on a particle in a discrete system can always be decomposed as a sum of 
central forces between particles, i.e., forces that are parallel to the lines connecting 
the particles. In other words, the strong law of action and reaction is always satisfied 
for such multi-body potentials. We show, that although the net force on a particle cal- 
culated using any extension is the same, its decomposition into central forces is gen- 
erally different for different extensions. Using this result we show that the pointwise 
stress tensor resulting from the Irving-Kirkwood-Noll procedure is non-unique and 
symmetric. We also show, that a generalization of Noll's lemmas [48] to non-straight 
bonds gives a non-symmetric stress tensor that may be important for particles with 
internal structure, such as liquid crystals. 

The macroscopic stress tensor corresponding to the pointwise stress tensor de- 
scribed above is obtained in Section [4] through a procedure of spatial averaging. The 
connection between this stress and the stress tensors obtained via the direct spatial 
averaging procedure introduced by Murdoch fl43jj47l and Hardy l23l is explored 
and in the process the Murdoch-Hardy procedure is systematized and generalized to 
multi-body potentials using the results of SectionfJ] The non-uniqueness of the stress 
tensor, inherent in the Murdoch-Hardy procedure is studied and a general class of 
possible definitions under this procedure are identified. The connection between the 
non-uniqueness in the Murdoch-Hardy procedure and the non-uniqueness mentioned 
in Section[3]is addressed. 

In Section|5] various stress definitions including the Hardy stress, the Tsai traction 
and the virial stress are shown to be special cases of the macroscopic stress tensor de- 
rived from the extended Irving-Kirkwood-Noll procedure in Section|4] The original 
definitions for these measures are generalized in this manner to multi-body potentials. 
The existence of different extensions for the potential energy function, which led to 
non-uniqueness of the pointwise stress tensor discussed in Section [3] also result in 
the non-uniqueness of these definitions. However it is shown that the difference in 
the macroscopic stress tensor resulting from this non-uniqueness tends to zero in the 
thermodynamic limi^ Another source of non-uniqueness explored in this section is 
that given any definition for the stress tensor, a new definition, which also satisfies the 
balance of linear momentum, can be obtained by adding to it an arbitrary tensor field 
with zero divergence. It is shown that in the thermodynamic limit the macroscopic 
stress tensor obtained in Section|4]converges to the virial stress derived in Section|2] 

To address practical aspects of the different definitions obtained within the unified 
framework, Section|6]describes several "numerical experiments" involving molecular 
dynamics and lattice statics. These simulations are designed to examine the behavior 
of these stress definitions, including their convergence with averaging domain size 
and their symmetry properties. Our conclusions and directions for future research are 
presented in Section|7] 



1 The thermodynamic limit is the state obtained as the number of particles, N, and the volume, V, of 
the system tend to infinity in such a way that the ratio N/V is constant. 
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Notation 

In this paper, vectors are denoted by lower case letters in bold font and tensors of 
higher order are denoted by capital letters in bold font. The tensor product of two 
vectors is denoted by the symbol "(g)" and the inner product of two vectors is denoted 
by a dot The inner product of two second-order tensors is denoted by A 
second-order tensor operating on a vector is denoted by juxtaposition, e.g., Tv. The 
gradient of a vector field, v(x), is denoted by V x v{x), which in indicial notation 
is given by [V x i;]y = dvi/dxj. The divergence of a tensor field, T(x), is denoted 
by div x T(x). The divergence of a vector field is defined as the trace of its gradient. 
The divergence of a second-order tensor field in indicial notation (with Einstein's 
summation convention) is given by [div^Tji = dTij/dxj. The notation described 
above is followed unless otherwise explicitly stated. 



2 Stress in an equilibrium system 

In this section, we obtain expressions for the Cauchy stress in an equilibrium sys- 
tem using the technique of canonical transformations. The basic philosophy behind 
canonical transformation is explained in the next section. 



2.1 Canonical transformations 

Consider a system consisting of N point masses whose behavior is governed by clas- 
sical mechanics. Let q a (t) and p a (t) (a = 1,2, ... , N) denote the generalized coor- 
dinates and momenta of the systemjj For brevity, we sometimes use q(t) andp(i) to 
denote the vectors (qi(t), q2(t), ■ ■ ■ , qjv(i)) and (pi(t),p2(t), . . . ,piy(t)), respec- 
tively. The time evolution of the system can be studied through three well-known 
approaches, referred to as the Newtonian formulation, the Lagrangian formulation, 
and the Hamiltonian formulation. The first approach is used in molecular dynamics 
simulations, while the latter two approaches are more elegant and can sometimes be 
used to obtain useful information from systems in the absence of closed-form solu- 
tions. 

In the Lagrangian formulation, a system is characterized by the vector q(t) and a 
Lagrangian function C, given by 

C{q,q;t)=T{q)-V{q), (2.1) 

where T is the kinetic energy of the system, V is the potential energy of the system, 
and q(t) represents the time derivative of q(t). It is useful to think of q as a point in 
a 3iV-dimensional configuration space. The time evolution of q(t) in configuration 
space is described by a variational principle called Hamilton's principle. Hamilton's 



2 In a general theory of canonical transformations, q a and p a need not denote the actual position and 
momentum of particle a. 
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principle states that the time evolution of q(t) corresponds to the extremum of the 
action integral defined as a functional of q by 

A[q] = ^ C(q,q;t)dt, (2.2) 

where t±, t^, q{t\) and qfa) are held fixed with respect to the class of variations 
being considered IT321 Section V.l]. In mathematical terms, we require that 

SA = 0, (2.3) 

while keeping the ends fixed as described above. The Euler-Lagrange equation aris- 
ing from ( 12.3b is 

d f dC\ dC 



- 0. (2.4) 

dt \oq a J oq a 

The Lagrangian formulation is commonly used as a calculation tool in solving simple 
problems. 

Next, we note that the Lagrangian is the Legendre transform of the Hamiltonian 
■H, ||32] Section VI.2], 

£{q,<i;t) = sup[p • q - H(p,q;t)]. (2.5) 
p 

The Hamiltonian is the total energy of the system. Using the Hamiltonian, equa- 
tion ( 12.3b can be rewritten as 

S I' \p-q-H(p,q;t)}dt = 0. (2.6) 
Jti 

Note that in ( 12.3b , the variation is only with respect to q, whereas in ( 12.6b , the func- 
tional depends on the functions q and p, and variations are taken with respect to 
both q and p independently. In both cases, t\, t%, q(t\) and qit-j) are held fixed. 
The variational principle given in ( 12.6b is commonly referred as the modified Hamil- 
ton's principle l22l or simply as the "Hamiltonian formulation". The advantage of 
the Hamiltonian formulation lies not in its use as a calculation tool, but rather in the 
deeper insight it affords into the formal structure of mechanics. The Euler-Lagrange 
equations associated with ( 12.6b are 

q a = y Pa n, (2.7) 
p a = -V qa H, (2.8) 

commonly called Hamilton's equations. The above equations are also referred to as 
the canonical equations of motion^ 

It is important to note that the Hamiltonian formulation is more general than the 
Lagrangian formulation, since it accords the coordinates and momenta independent 
status, thus providing the analyst with far greater freedom in selecting generalized 



3 The term "canonical" in this context has nothing to do with the canonical ensemble of statistical 
mechanics. The terminology was introduced by Jacobi to indicate that Hamilton's equations constitute the 
simplest form of the equations of motion. 
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coordinates. We now think of (q, p) as a point in a 6iV-dimensional phase space, as 
opposed to the 3iV-dimensional configuration space of the Lagrangian formulation. 
The choice of q and p is not arbitrary, however, since the selected variables must 
satisfy the canonical equations of motion. For this reason q and p are called canonical 
variables. 

The requirement that the generalized coordinates and momenta must be canoni- 
cal means that new sets of generalized coordinates can be derived from a given set 
through a special kind of transformation defined below. 

Definition 1 Any transformation of generalized coordinates that preserves the canon- 
ical form of Hamilton's equations is said to be a canonical transformation^ 

The construction of canonical transformations is facilitated by the introduction of 
generating functions as explained below. 

Generating functions 

Consider two sets of canonical variables (Q,P) and (q,p), related to each other 
through a canonical transformation given by 

Q = Q(q,P,t), P = P(q,p,t). (2.9) 
Since the variables are canonical, they satisfy the modified Hamilton's principle in 
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6 [ 2 [p-q-H(p,q;t)]dt = Q, (2.10) 

[ 2 \p-Q-U(P,Q;t) 

Ju 1 J 



r«2 r 

' dt = 0, (2.11) 



where H is defined later as part of the canonical transformation. The integrands of 
( 12.10b and ( 12.111 ) can therefore only differ by a quantity whose variation after inte- 
gration is identically zero. A possible solution is 

f' 2 r i f t2 dC 

5 p-q-P-Q-(H-H) dt = S —dt, (2.12) 
Jti L J Jti dt 

where G is an arbitrary scalar function of the canonical variables and time, with 
continuous second derivatives. The integral on the right is only evaluated at fixed 
integration bounds and its variation is zero. This is not obvious since there is no 
restriction on the variation of the momenta at the ends. We assume this to be true to 
avoid the introduction of differential forms. For a mathematically rigorous argument 
refer to fT\ Section 450. The difference between the integrands of ( I2.10l i and ( 12.1 II ) 
therefore satisfies, 

dG-p-dq + P-dQ + (H-n)dt = 0. (2.13) 



4 This definition suffices for our purpose, but a more correct definition can be found in (TJ using 
differential forms. 

5 Briefly the proof is based on the symmetry present in the geometry of any Hamiltonian system 
commonly called symplectic geometry. 
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Now, consider the case where G = Gi(q, Q, t). The total differential of G is then 

d_G_ 

Bt 



BG 

dG = \7 q G 1 -dq + \7 Q G 1 -dQ + —±dt. (2.14) 



Substituting ( 12.14b into ( 12.131 ) gives 

(VqGi - p) ■ dq + (V Q Gi +P)-dQ+[—^-+n-n)dt = Q. (2.15) 



Since q, Q and t are independent, the above equation is satisfied provided that 

P a = 1 r ± , P a = -^pr, H=H + —±. (2.16) 
Bq a oQ a at 

The above relations define the canonical transformation. Since G\ generates the 
transformation, it is commonly called the generating function of the canonical trans- 
formation. Note that if G\ does not depend on time t, then H = Ti. 

The generating functions of the form, G = Gi(q,Q,t), does not generate all pos- 
sible canonical transformations. In general, there are four primary classes of generat- 
ing functions where the functional dependence is (q, Q), (q, P), (p, Q) and (p, P)@ 
We have already encountered the first class, where G = Gi(q, Q, t). The remaining 
classes can be obtained from the first through Legendre transformations. Consider for 
example, the following definition, 

G = G 3 (p,Q,t)+q-p. (2.17) 

The total differential of this expression is 

BG 

dG = V„G 3 • dp + V Q G 3 • dQ + -^-dt + q ■ dp + p ■ dq. (2.18) 

Bt 

Substituting the above equation into (12.131 ) gives 

(V P G 3 + q) ■ dp + (V Q G 3 + P) ■ dQ + + H - Hj dt = 0, (2.19) 

which leads to the following canonical transformation: 

Qa = -f^, Pa = -^-, n=H + ^. (2.20) 

Bp a BQ a Bt 

The other two classes of transformation can be derived in a similar way. 

Finally, an important property of a canonical transformation is that it preserves 
the volume of any element in phase space, i.e., dqdp = dQdP E2l page 402]. This 
means that for a change of variables between (p, q) and (P, Q), the Jacobian of the 
transformation is unity. 



6 In addition to these four classes of transformation, it is possible to have a mixed dependence, where 
each degree of freedom can belong to a different class 1221 . 
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2.2 A derivation of the stress tensor under equilibrium conditions 

In this section, we use the method of canonical transformations to derive an expres- 
sion for the Cauchy stress tensor. In continuum mechanics, a body is identified with 
a regular region of Euclidean space £ referred to as the reference configuration. Any 
point X G B is referred to as a material point. The body B is deformed via a smooth, 
one-to-one mapping tp : £ —> £, which maps each X G B to a point, 

x = tp{X), (2.21) 

in the deformed configuration^ where we have assumed that the deformation is inde- 
pendent of time. The deformation gradient F is defined as 

F(X) = V X ifi. (2.22) 

The mapping ip is assumed to satisfy the condition that dot F is strictly positive. The 
Cauchy stress, <x, is defined by 11371 

<r(T,F) = -^V F il>F T , (2.23) 

where ip(T, F) is the Helmholtz free energy density function relative to the reference 
configuration. We are only focusing on a conservative elastic body. 

A system in thermodynamic equilibrium^ can by definition only support a uni- 
form state of deformation. Therefore, our material system is deformed via the affine 
mapping 

q a = FQ a . (2.24) 



7 We adopt the continuum mechanics convention of denoting variables in the reference configuration 
with upper-case letters, and variables in the deformed configuration with lower-case letters. 

8 A system is said to be in a state of thermodynamic equilibrium when all of its properties are inde- 
pendent of time and all of its intensive properties are independent of position 1651 . To stress this, the term 
uniform state of thermodynamic equilibrium is sometimes used to describe this state. 

9 To understand this mapping, consider a system of N particles with positions q a (a = 1, 2, . . . , N) 
confined to a parallelepiped container defined by the three linearly independent vectors l\, I2 and I3, 
which need not be orthogonal. This selection is done for convenience and does not limit the generality of 
the derivation as explained below. The position of a particle in the container can be expressed in terms of 
scaled coordinates £ [0, 1] as 

q«=&k> (*) 
where Einstein's summation convention is applied to spatial indices. The deformation of the container is 
defined relative to a reference configuration where the cell vectors are L 1 , L2 and L;$ . The current and 
reference cell vectors are related through an affine mapping defined by F, 

li = FLi. (**) 

Equations Q and (**\ can be combined to relate the position q a of particle a in the deformed configura- 
tion with its position in the reference configuration Q a , 

q a =tf(FL i ) = F(gL i ) = FQ a . (***) 



This is exactly the mapping defined in )2.24t . It provides a direct relationship between the positions of 
particles in the reference configuration and their position in the deformed configuration. Note that the 
assumed (parallelepiped) shape of the container does not enter into the relation, q a = FQ a , which 
means that this relation holds for a container of any shape. 
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It is clear that if we enforce this mapping on our system, with no change in the 
momentum coordinates, then the newly obtained variables will not satisfy Hamil- 
ton's equations. Therefore any change of variables should be governed by a canoni- 
cal transformation. The following generator function provides the desired canonical 
transformation 

G 3 (p,Q) = -^2p a -FQ a . (2.25) 

a 

Substituting this generating function into ( 12.20b gives 

q a = -^ = FQ a , Pa = -^ = F T p a , n = H. (2.26) 

The first relation in the above equation is the desired transformation in (12.24-b . The 
second relation is the corresponding transformation that the momentum degrees of 
freedom must satisfy, so that the new set of coordinates (Q, P) are canonical. The 
third relation refers to the Hamiltonian of the system, which is assumed to be given 
by 

N 

H(p,q) = Y, E Z J ^ + V(<lu..-,qN), (2.27) 

a— 1 

where V denotes the potential energy of the system. Expressed in terms of the refer- 
ence variables, (12.27b becomes 



H(P, Q, F) = H(p(Q, P, F),q(Q, P, Fj) 
F T P a F T P C 



E 



2m c 

a— 1 



V(FQ 1 ,...,FQ N ). (2.28) 



We now proceed to derive the expression for the Cauchy stress tensor using ( 12.231 ). 
The Helmholtz free energy density for the canonical ensemble is given by (26) 



^T,F) = -M^, (2.29) 
vo 

where ks is the Boltzmann's constant, T is the absolute temperature, Vq is the vol- 
ume of the body in the reference configuration, and Z(T, F) is the partition function 
defined as 

Z{T, F) := J e-^ k - T dPdQ, (2.30) 



It is important to note that does not impose a kinematic constraint that dictates the position of 

particle a in the deformed configuration based on its position in the reference configuration (as does the 
Cauchy-Born rale used in multiscale methods 1581 ). We will see later that this will merely be used as a 
change of variables, where, instead of integrating over the deformed configuration with the variables q, 
the integration is carried out over a given reference configuration using the variables Q. In both cases the 
same result is obtained. However, by using the referential variables the dependence on the deformation 
gradient is made explicit. 
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where h is Planck's constant and Jo denotes the phase space in the reference con- 
figuration. With this definition, the statistical mechanics phase average of a function 
A(P, Q) in the canonical ensemble is 

(A)(T,F)= f A(P,Q)W c (P,Q,T,F)dPdQ, (2.31) 

Jr n 



where 



W C (P,Q,T,F) 



1 



,-H(P,Q,F)/k B T 



(2.32) 



N\h 3N Z 

is the canonical distribution function. Substituting ( 12.291 ) and ( 12.301 ) into ( 12.231 ). we 
obtain 



k B T 



; v f zf t = - h F n)F T , 



(det F)V a Z"' V 
where in the last step we have used the identity V = (det F)Vq and where 



(2.33) 



1 



N\h 3N 
1 



k R TN!h 3N 



e -n/k B T d Q dp 
V F 'He-' fL/kBT dQdP. 



(2.34) 



Next, we compute V f%. In our derivation, we make use of indicial notation and 
the Einstein summation rule. To accommodate for the spatial indices, we push a 
representing the particle to the superscript position. Following this adjustment, we 
have 







&H 



dF u dFt 



u 



E 



2m a 



V(q\...,q N ) 



E 



m a dFu k 



From ( 12.26b . we have 



dV dg% ' 
dq% dF u 

(2.35) 



dq% 


d 


dFu 


dFu 


dvt 


d 


dFu 


dFu 



(F kL Q2) = S lk Q 



( F Lk P L) 



J ! 



r Jk r Li r L 



where in ( 12.37b . we have used the following identity: 



dF 



U 



p- 1 p- 1 
~ r Li r Jk ■ 



Substituting (12.36b and (12.37b into ( 12.35b . we have 



r Jk Pt ' 



(2.36) 
(2.37) 

(2.38) 



&H 

dF, 



u 



E 



\P?FJkPt 



Ja,iWj 



(2.39) 



13 



where = —dV/dr a is the internal force, defined in the deformed configuration, 
on particle a In direct notation, we have 



p a ®F 1 p a , 



fa ® Q<* 



(2.40) 



Substituting the above equation into ( I2.33l l and using (12.26b , we obtain an expression 
for the Cauchy stress: 

«x(T, F) = - I V / + ® g a \ , (2.41) 



where the phase averaging is now being performed with respect to the variables p and 
q. The switch from phase averaging over P and Q in d2.311 l to p and q above can be 
made because canonical transformations preserve the volume element in phase space 
as explained at the end of Section lXTl 

The expression in ( 12.411 ) for the Cauchy stress tensor is called the virial stress. A 
simpler derivation of the virial stress, based on time averages, is given in Appendix 
lAl Although, the derivation here made use of the canonical ensemble, it is expected to 
apply to any ensemble in the thermodynamic limit (see footnote[TJon page|5j where 
all ensembles are equivalent. Continuum mechanics also tells us that the Cauchy 
stress tensor is symmetric, something that is not evident from the above equation. 
Discussion of the symmetry of the stress tensor, which hinges on an important prop- 
erty of f™ 1 , is postponed to Section|5] 

The virial stress defined above corresponds to the macroscopic stress tensor only 
under conditions of thermodynamic equilibrium in the thermodynamic limit. We now 
show that this expression for the stress tensor, as well as all other expressions in 
common use, can be derived as limiting cases of a more general formulation which 
begins with the Irving-Kirkwood-Noll procedure. We refer to this as the "unified 
framework" for the stress tensor. 



3 Continuum fields as phase averages 

In this section, we discuss the Irving and Kirkwood procedure (27), which laid the 
foundation for the microscopic definition of continuum fields for non-equilibrium 
systems. This work was later extended by Walter Noll [480, who showed how 
closed-form analytical solutions can be obtained for the definition of certain contin- 
uum fields, which otherwise involved a non-rigorou^3 series expansion of the Dirac 

10 There is a subtle point here. Since we are using the canonical ensemble, the Hamiltonian H neglects 
the interaction term of the system with the surrounding "heat bath". This means that the potential energy 
V in H only includes the internal energy of the system and, therefore, its derivative with respect to the 
position of particle a gives the force /^ nt on this particle due to its interactions with other particles in the 
system. 

"An English translation of this article appears in the current issue of the Journal of Elasticity. 

12 The derivation is non-rigorous in the sense that expressing the stress tensor as a series expansion is 
only possible when the probability density function, which is used in the derivation, is an analytic function 
of the spatial variables 1481 . 
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delta distribution in the original procedure. We refer to the procedure proposed by 
Noll in [48 1 as the Irving-Kirkwood-Noll procedure. The derivation presented in this 
section largely follows that of Noll l48l . but extends it to more general atomistic 
models. 

Consider a system Ai modeled as a collection of N point masses/particles, each 
particle referred to as a (a = 1,2, ... , N). We use the terms "particle" and "atom" 
interchangeably. The position, mass and velocity of particle a are denoted by x a , 
m a and v a , respectively. The complete microscopic state of the system is known, at 
any instant of time, from the knowledge of position and velocity of each particle in 
R 3 . Hence, the state of the system at time t, may be represented by a point E(t) in 
a 6iV-dimensional phase space0. Let r denote the phase space. Therefore any point 
S{t) € r, can be represented as, 



In reality, the microscopic state of the system is never known to us, and the only 
observables identified are the macroscopic fields as defined in continuum mechanics. 
We identify the continuum fields with macroscopic observables obtained in a two- 
step process: (1) apointwise field is obtained as a statistical mechanics phase average; 
(2) a macroscopic field is obtained as a spatial average over the pointwise field. The 
phase averaging in step (1) is done with respect to a probability density function W : 
r x M + -> M + of class C 1 defined on all phase space for all t (W c , defined in ( 12.32b . 
is an example of a stationary (time-independent) probability density function defined 
for the canonical ensemble). The explicit dependence of W on time t, indicates that 
our system need not be in thermodynamic equilibrium. 

As discussed in Section [2] the evolution of S(t) in the phase space is given by 
the following set of 2N first-order equations (Hamilton's equations of ( I2.7b - (l2.81 l): 



where p := (pi,j>2, ■ • ■ ,Pn), Pa denotes the momentum of each particle, and 
H(p, x) is the Hamiltonian of the system. 

The basic idea behind the original Irving and Kirkwood procedure is to pre- 
scribe/derive microscopic definitions for continuum fields, such that they are con- 
sistent with the balance laws of mass, momentum and energy. To arrive at these defi- 
nitions, we repeatedly use the following theorem, commonly referred to as Liouville 's 
theorem, which relates to the conservation of volume in phase space. 

As a system evolves, the phase space _T is mapped into itself at every instant 
of time, and this mapping is governed by ( 13.21 i. If g t denotes this mapping, then 
Liouville's theorem essentially says that for any subset U of r, the volume of U 
remains invariant under the mapping g t . This can be be formally stated as, 

13 The usual convention is to represent the phase space via positions and momenta of the particles. For 
convenience, in this section, we represent the phase space via positions and velocities of the particles. 



3(t) = (x 1 (t),x 2 (t),. . . ,x N (t);v 1 (t),v 2 (t), . . -,v N (t)) 

=■■ (*(*);«(*))■ 



(3.1) 



P = -v x h 
x = V P H, 



(3.2a) 
(3.2b) 
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Liouville's Theorem For any t/Cf, volume is preserved under the one-parameter 
group of transformations of phase space, g t : U — > r, given by the mapping 

(as(0),p(0)) ^ (x(t),p(t)), 
where x(t) and p(t) are solutions of the Hamilton's system of equations <\3.2\ , i.e., 

vol(U)=vol(g t U). (3.3) 



Proof Let vol(g t U) denote the material time derivative of vo\(g t U) in the sense that 
3(0) is held fixed while performing this differentiation. Then we have, 



vo\{g t U)= / d3(t) = / (detF)dS , 

JgtU JU 

where F(3 ,t) := V Sg S(S ,t), 3(3 ,t) = g t (B ) and S = 3(0). Using the 
fact that 



dctF = (dctF)tr(FF- 1 ), 



we obtain 



Let 



vol(g t U)= / (detF)ti(FF- 1 )dS„. 
Ju 



From chain rule, we have 



±_ d{VS) 



dt 



=o=g t l (3) 



V S S = FF 



5 =g t 

Therefore div 3 = tr(FF~ 1 ). Equation ( 13.4b can now be rewritten as 



vo\(g t U) = / (detF)(div 3) \ s{t)=gt(So ) d3 Q . 
Ju 

But from ( 13.11 ) and ( 13.21 i we also have, 

div 3 = div x x + div p p = div^(V p H) — div p (V x 'H) = 0. 



(3.4) 



(3.5) 



(3.6) 



(3.7) 



Therefore vo\(g t U) = for arbitrary t. Thus ( 13.3l l holds. 



have 



Let W(3; t) denote the probability density function defined on gt{r). Hence, we 



/ W(S(t);t)dS(t) = f W(3 Q ;0)(detF)dS . 

JgtU JU 



As a consequence of Liouville's theorem, we have, det F = 1. Therefore 

d 



4 / W(S(t);t)dS(t)=0. 

dt JgtU 



(3.8) 



(3.9) 
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Since ( 13.91 ) holds for all U C r, we have VT(S(i); i) = 0. Hence, the time evolution 
of the probability density function is given by 



dW 



N 



tj VV X — > 

— + K- v,^ + v a ■ V Va W] = 0. 



(3.10) 



The above equation can be rewritten as 



dW " 



v x v 



0. 



(3.11) 



where, as before, V(x\,X2 1 ■ ■ ■ ,Xn) denotes the potential energy of the system. 
Equation ( 13.1 U is called Liouville's equation. 



3.1 Phase averaging 

Under the Irving-Kirkwood-Noll procedure, pointwise fields are defined as phase 
averages. This phase averaging is expressed via weighted marginal densities. For 
example, the pointwise mass density field is defined as 

p(x,t):=y m a / WS(x a — x) dxdv, (3.12) 

„ Jm 3N xm? N 

where the integral represents a marginal density defined on R 3 , S denotes the Dirac 
delta distribution, and V denotes summation from a = 1 to N. To avoid the Dirac 
delta distribution and for greater clarity we adopt Noll's notation as originally used 
in ll48l . Hence ( 13.121 i can be rewritten as 

P (x,t) = j:m a fwdx 1 ...dx a _ ld x a+1 ...dx Nd v 



E 



a (W\x a = x), (3.13) 



where (W \ x a = x) denotes an integral of W over all its arguments except x a and 
x a is substituted with x. Now consider the continuum velocity field. Unlike the defi- 
nition of pointwise density field, which appears unambiguous, the pointwise velocity 
field can be defined in different ways. It may seem more natural to define the contin- 
uum velocity in an analogous fashion to the density field, i.e., 

Ea (W\x a = X) 

Alternatively, the pointwise velocity field can be defined via the momentum density 
field, p(x, t), as follows: 

p(x, t) := 2^ m a (Wv a I x a = x) , (3.15) 

a 

v(x,t):= (3.16) 
p(x,t) 
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Note that definitions ( 13.14l i and ( 13.161 ) are equivalent for a single species material, 
but are not so in general. The definition given by ( 13.16b is the one used in practice. 
There are two reasons for this. First, the definition in ( 13.161 ) makes more physical 
sense since, following spatial averaging, it associates the continuum velocity with the 
velocity of the center of mass of the system of particles. Second, the definition in 
d3.161l satisfies the continuity equation as shown in Section [331 whereas ( 13.141 ) does 
notQ 



3.2 Regularity assumptions for the probability density function 

It is clear from the definitions in ( 13.13b , ( 13.15b and (13.161 ) that the integrals in these 
equations converge under appropriate decay conditions on W. The following two 
conditions are sufficient for the convergence of all the integrals and the validity of the 
results in this section 11481 : 

1 . There exists a 6 > such that the function 

JV N 

W(S;t)l[\\x a \\ 3 + s l[\\v p \\ 6 + s (3.17) 

o=l 0=1 

and its first derivatives are bounded by a constant that only depends on time. 

2. V(xi, x<2, ■ ■ ■ xjv) is abounded C 1 function defined on the phase space, and hav- 
ing bounded first derivatives^ 

Conditions (Q~|) and (fjji ensure the convergence of all the integrals considered in this 
section and swapping of integration and differentiation. Furthermore, let G(S; t) be 
any vector or tensor-valued function of class C 1 defined on the phase space for all t, 
and which, for suitable functions g(t) and h(t), satisfies the condition 

N 

sup (||G||, || div Ua G||, || div XQ G||) < g(t) TT \\v p \\ 3 + h(t), 

(3.18) 

where || ■ || refers to the norm defined through the inner product. Since the space of 
all tensors has a natural inner product defined as 

S :T = tr(S T T), (3.19) 



14 It would be interesting to explore how the equation of continuity fails for the definition in )3.14t 
by identifying the regions that act as sinks and sources. This is difficult to do for a general N particle 
system because the continuity equation quickly becomes unwieldy. Even for the much simpler case of a 
two-particle system, the answer is not trivial. A quick examination shows that the distribution of sinks and 
sources depends not only on the ratio of the masses but also on the probability density function W. 

15 If any two particles overlap, we would normally expect V — > oo. By specifying additional decay 
conditions for W, the case of unbounded V can be handled. For simplicity, we assume V to be bounded. 
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we have || S\\ = V S : S. Under these conditions on G(S; t), we havJ^l 

/ G-\7 Xa Wdx a = - [ Wdiv Xa Gdx a , (3.20a) 

JR3 JR3 

/ G-\7 Va Wdv a = - I Wdxv Va Gdv a . (3.20b) 

JR3 JR3 

The above identities are repeatedly used in deriving the equation of continuity and 
the equation of motion in the following sections. 



3.3 Equation of continuity 



Let us demonstrate that the pointwise fields defined in Section BTI satisfv the equation 
of continuity. The equation of continuity from continuum mechanics is given by Il37ll 



dp 

at 



+ div x (pv) = 0. 



(3.21) 



dp, x / dW 

(x,t) = ^m a { — 



From ( 13.131 ) we have 

a ' 

Using Liouville's equation in (13. 1 lb . we have 

a \ ^ 



m/3 



Now, consider the summand on the right-hand side of the above equation for a fixed 

a. From O.20bt . it is clear that ^ • V Vp W x a = x) = 0, for /3 = 1, 2, • • • N, 

and from ( I3.20al i. we also have (yp ■ V Xfi W \ x a = x) = 0, for j3 ^ a. Therefore 
the above equation simplifies to 



dp 
lit 



(x, t) = - rn Q (v Q • V Xq iy | a; Q 



Using the identity, 



div x (aw) = V^a • w, 



(3.22) 



where a(:c) is any C 1 scalar function of x, and u; is any vector independent of x, we 
obtain 

*^{x,t) = ~y] m a div x (VKu Q | sCq, = x) . 



16 If G is a second-order tensor or higher, then the dot product indicates tensor operating on a vector. 
Note that in 13. 201 . in the interest of brevity, we are breaking our notation of denoting a second-order 
tensor operating on a vector by juxtaposition. 
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Using ( I3.15l l and ( 13.161 ) for the definition of the pointwise momentum density field, 
we have 

^(x,t)+div x (pv) = 0, 

which is the continuity equation. We have established that the definitions given in 
( 13.12b and ( 13.15b identically satisfy conservation of mass. 



3.4 Equation of Motion 

The equation of motion from continuum mechanics is given by (37 
d(pv) 



dt 



d'w x {pv (g) v) = div x cr + b. 



(3.23) 



Here we identify cr with the pointwise stress tensor. From ( 13.151 ). we have 

dp, s I dW 



Again, using ( 13.111 ) we obtain, 
dp 



(x, t) = ^ m a ( V a ^ I -V Xfl W ■ Vp 



TO/3 



v,„v 



TO/3 



(3.24) 



Now, consider the summand on the right-hand side of the above equation for fixed a 
and j3. Using d3.20ab . we have ((v a ® vp)\7 Xl) W \ x a = x) = 0, for (3 ^ a. From 
(13.20M . we have ((v a ® V Xli V)V V/3 W \ x a = x) = 0, for (3 ^ a, and for /3 = a, 
we have 

(( Va ® V Xo V)V„ Q W | SB a = x) = - (V x „ VW^ | 35a = SC) , 

using the fact that div u (it ® w) = w, for any vector u and for any vector w inde- 
pendent of it. Therefore (13.241 ) simplifies to 

T^OM) = ~X] ma! (( Ua! ^cOVa^W I =X) - ^2(WV Xa V \x a =x). 

a a 

(3.25) 

Using the identity, 

div x (aT) = TV x a, (3.26) 

where a(x) is any C 1 scalar function of x, and T is any tensor independent of x, we 
can rewrite (13.25b as 

^■{x,t) = - div x ^2m a ((v a ® v a )W \ x a = x) - ^ (WVaj a V | aj„ = as) . 

a. a 

(3.27) 
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Now, note that the term v a ® v a can be written as 

v a ®v a = (v a — v) ® («„ -i7)+i7(g>i7 Q .-|-i7 Q: (gii7 — I7g)i7 

= < o1 (8) < cl + v <g> 17 Q + u Q ® u - u <g> 17, (3.28) 

where i7jj eI is the velocity of particle a relative to the pointwise velocity field. Con- 
sider the first term on the right-hand side of ( 13.27b . Substituting (13.28b into this ex- 
pression we have, 

- div x 2J m a ((«£» ® UaJW | SC a = x) 

a 

= - ^ m Q div x ((17" 1 ® i7jj Cl )W | x a = sc) - div x 2J [17 <g) m a (i> Q VK | a5 a = x) 

a a 

+ m a (v a W | x a = x) (g) 17 — TOq, (W \x a = x)v®v] 
= - div x ^ m a («' <g> <')W | x a = x) - div x (pi; <g> 17), (3.29) 

a 

where we have used ( 13.13b . ( 13.15b and ( 13.16b in the last step. Substituting ( 13. 29b into 
(13.27b . we obtain 

*^-{x,t) +div x (pv ®v) = -^m a div x ((v r f®v T ( ?)W \ x a = x) 

a. 

-^2(WV Wet V\x a = x). (3-30) 

ot 

The left-hand sides of (13.30b and (13.23b are identical. Therefore, the right-hand sides 
must also be equal. Hence 

div x a+b = - J2 m a div x <« cl ® vf)W | x a = x)-J2 {WV Xa V \x a = x). 

a a 

(3.31) 

To proceed, we divide the potential energy V(xi, X2, ■ ■ ■ , xn) into two parts: 

1. An external part, V cx t> associated with long-range interactions such as gravity or 
electromagnetic fields. 

2. An internal part, V; n t, associated with short-range particle interactions. 

It is natural to associate V cx t with the body force field b in (13.31b . We therefore define 

b(x, t) as 

b(x, t) := - (WV Xa V cxt | x a = x) . (3.32) 

a 

Substituting (13.32b into ( 13.31b . we have 

div x <t = - ^2 ™« div * (« el ® v T f)W | x a = x)-^2{WV Xa V iv .t \x a =x). 

a a 

(3.33) 
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From ( 13.33b . we see that the pointwise stress tensor has two contributions: 

tr{x,t) = <Tk{x,t) + <r v (x,t), (3.34) 

where cr^ and <x v are, respectively, the kinetic and potential parts of the pointwise 
stress. The kinetic part is given by 

<r k (x,t) = -J2 m a < « o1 <g> vf)W \x a = x). (3.35) 

a 

It is evident that the kinetic part of the stress tensor is symmetric. The presence of 
a kinetic contribution to the stress tensor appears at odds with the continuum defi- 
nition of stress that is stated solely in terms of the forces acting between different 
parts of the body. This discrepancy has led to controversy in the past about whether 
the kinetic term belongs in the stress definition ll67ll . The confusion is related to the 
difference between absolute velocity and relative velocity defined in ( 13.28b Il63l . The 
kinetic stress reflects the momentum flux associated with the vibrational kinetic en- 
ergy portion of the internal energy. 

Continuing with ( 13.33b . the potential part of the stress must satisfy the following 
differential equation: 

div x <r v (x,t) = { w f'T \x a =x), (3.36) 

a 

where 

f'T := -V^Vmt, (3.37) 

is the force on particle a due to internal interactions. Equation ( 13.361 ) needs to be 
solved in order to obtain an explicit form for cr v . In the original paper of Irving and 
Kirkwood |271, this was done by applying a Taylor expansion to the Dirac delta dis- 
tribution appearing in the right-hand side of the equation. In contrast, Noll showed 
that a closed-form solution for cr v can be obtained by recasting the right-hand side 
in a different form and applying a lemma proved in |48l . We proceed with Noll's 
approach, except we place no restriction on the nature of the interatomic potential 
energy Vint- The potential energy considered in l27l and l48l is limited to pair poten- 
tials. 

General interatomic potentials 

In general, the internal part of the potential energy, also called the interatomic poten- 
tial energy, depends on the positions of all particles in the system: 

Vua=V} J a(x 1 ,X2,...,xn), (3.38) 

where the "hat" indicates that the functional dependence is on absolute particle posi- 
tions (as opposed to distances later on). We assume that Vint : K 3Ar — > M is a contin- 
uously differentiable function[3 This function must satisfy the following invariance 
principle: 



17 Note that this assumption may fail in systems undergoing first-order magnetic or electronic phase 
transformations. 
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The internal energy of a material system is invariant with respect to the Eu- 
clidean group Q ;= {x i-» Qx + c \ x G R 3 , Q € 0(3), c e K 3 }, where 
0(3) denotes the full orthogonal group. 

To exploit this invariance, let us consider the action of Q on R 3JV , i.e., the action of any 
combination of translation and rotation (proper or improper), which is represented by 
an element g : x M> Qx + c in Q, on any configuration of N particles represented by 
a vector (x t , . . . ,x N ) e R 3W : 

.g • (x 1} x N ) = (Qxi + c, . . . , Qx N + c). (3.39) 

This action splits R 3W into disjoint sets of equivalence classes IPT41 . which we now 
describe. For any u = (x-i , . . . , x^) G R 3Ar , let O u C M. 3N denote an equivalence 
class which is defined asF^I 

O u :={g-u\geg}, (3.40) 

where g ■ u denotes the action of g on u defined in ( 13.39b . In other words, O u repre- 
sents the set of all configurations which are related to the configuration u by a rigid 
body motion and/or reflection. Due to the invariance of the potential energy, we can 
view the function Vint as a function on the set of equivalence classes, i.e., 

V int (O u ) = V int (w), (3.41) 

because 

V in t («) -V int («) VveO u . (3.42) 

Now, consider a set S C R^f^ -1 )/ 2 , defined as 

S := {(ri2,ri3,...,rijv,r 2 3, ...,r(jv_i)jv) I 

r Q/3 = \\x a -x f3 \\,(x 1 ,...,x N ) e R 3JV }. (3.43) 

In other words, the set 5 consists of all possible N(N — l)/2-tuples of real numbers 
which correspond to the distances between N particles in R 3 ^J In technical terms, 
the coordinates of any point in S are said to be embeddable in R 3 . Note that S is a 
proper subset of K JV ( JV - 1 )/ 2 as it consists of only those N(N — l)/2-tuple distances 
which satisfy certain geometric constraints. In fact, the set S represents a (37V — 6)- 
dimensional manifold in K™!^ 1 )/^ commonly referred to as the shape space. 

Let 4> be the mapping taking a point in configuration space to the corresponding 
set of distances in S, i.e., : R 3JV — s- S : (xx, ■ ■ ■ ,a;jv) H> (^12, • • • , ^(iv-i)jv). 
where r a p from here onwards is used to denote \\x a — Xp\\, Since the Euclidean 
group preserves distances, it immediately follows that the map 

4> : {Equivalence classes} — > S, (3.44) 

18 The notation "{g ■ u \ g 6 Q}" should be read as "the set of all g ■ u, such that g is in the Euclidean 
group Q". 

"The key here is that not all N(N — l)/2 combinations of real numbers constitute a valid set of physi- 
cal distances. The distances must satisfy certain geometric constraints in order to be physically meaningful 
as explained below. 
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defined as <f>(O u ) = <j>{u), is a bijection (one-to-one and onto mapping) from the 
set of equivalence classes to the set This essentially means that for every set of 
equivalent configurations, i.e., configurations related to each other by a rigid body 
motion and/or reflection, there exists a unique N(N — l)/2-tuple of distances and 
vice versa. From ( 13.4 U and ( 13.441 ). it immediately follows that the potential energy 
of the system can be completely described by a function Vint : 5 — > M, defined as 

Vint(s) := Vtat(0 -1 (*)) Vs G S. (3.45) 



We now restrict our discussion to those systems for which there exists a contin- 
uously differentiable extension of V nt , defined on the shape space, to ]R Ar ( 7V_1 )/ 2 |^] 
This is justifiable because of the fact that all interatomic potentials used in prac- 
tice, for a system of N particles, are either continuously differentiable functions on 
jn>iv(iV-i)/2^ or can eas iiy be extended to one. For example, the pair potential and the 
embedded-atom method (EAM) potential |[T2l are continuously differentiable func- 
tions on R N ( N ~ 1 )/ 2 , while the Stillinger- Weber El and the Tersoff J60) potentials 
can be easily extended to R^f"- 1 )/ 2 by expressing the angles appearing in them 
as a function of distances between particles. Therefore, we assume that there exist a 
continuously differentiable function Vint : R^f^" 1 )/ 2 —> ]R, such that the restriction 
of Vmt to S is equal to V; n t: 

V nt 0) = Vint («) Vs = (ri 2 ,...,r (JV -i)jv) G S. (3.46) 

An immediate question that arises is whether this extension is unique in a neigh- 
borhood of s G S. Note that for TV < 4, 3N — 6 = N(N - l)/2. Therefore, for 
N < 4, for every point seS, there exists a neig hborhood in M.n(n-i)/2 which Hes 
in S. However, for N > 4, there may be multiple extensions of Vi n t. 

As noted above, the reason we are considering an extension is to define the par- 
tial derivative of the potential energy with respect to each coordinate of a point in 
jjAf(jv-i)/2 -phis will be used later to define the stress tensor. For example, the 
partial derivative of Vi n t(Ci2, • ■ • , C(n-i)n) with respect to £12 at any point s = 
{ru, r(jv_i)jv) G S, defined as 

dVint i , _ y Vint(ri2 + £, ■ ■ ■ ,T-]y(jV-l)/2) ~ Vjnt(?"l2, ■ ■ ■ , ? N -\) j 2) 



(3.47) 



requires us to evaluate the function at non-embeddable points. 



4> is surjective (onto) by the definition of S. The proof that it is injective (one-to-one) is similar to 
the proof of the basic imariance theorem for the simultaneous invariants of vectors due to Cauchy, which 
can be found in (62) Section 11]. 

21 The extension is necessary since Vi nt is defined in )3.45t only on the set S. We need to extend the 
definition to all points in R^^ -1 )/ 2 , whether they correspond to a set of physical distances or not, in 
order to be able to compute derivatives as explained later in the text. This issue has been overlooked in 
the past (see for example |13| ). which leads to the conclusion that the stress tensor is always symmetric. 
It turns out that this conclusion is correct (at least for point masses without internal structure), but the 
reasoning is more involved as we show later. 



24 



It will be shown later that the quantity evaluated in (13.47b may differ for different 
extensions. On the other hand, V Xa V[ n t is uniquely defined for any extension. This 
is because 

= V a!o Vint(^" 1 (s)) 

= V. a Vint(ti), (3.48) 

where = O u , which implies <p(u) = s@ and we have used (13.46b . (13.45b 

and ( 13.41b in the first, second and the last equality respectively. 

We next address the possibility of having multiple extensions for the potential en- 
ergy by studying the various constraints that the distances between particles have to 
satisfy in order to be embeddable in R 3 . We demonstrate, through a simple example, 
how multiple extensions for the potential energy can lead to a non-unique decompo- 
sition of the force on a particle, which in turn leads to a non-unique pointwise stress 
tensor. 



Central-force decomposition and the possibility of alternate extensions 

We will now show that the force on a particle can always be decomposed as a sum of 
central forces. The force on a particle due to internal interactions is defined in ( 13.37b . 
Therefore, for any configuration u £ M. 3N , we have 

/i nt H = -V Xa V int (u). (3.49) 

Using ( 13.481 ). ( 13.491 ) takes the form 

= ^ f a p(u), (3.50) 

13 

where s = <j>{u) = (r 12 , . • . , r (N _ 1)N ) and 



dVu 



■(0(«))^^ifa</3, 



^^ : -{g(^))4^ifa> A (151) 

is the contribution to the force on particle a due to the presence of particle (3. 

Note that f a p is parallel to the direction x@ — x a and satisfies f a p = —fp a - 
We therefore note the important result that the internal force on a particle, for any 
interatomic potential that has a continuously differentiable extension, can always be 
decomposed as a sum of central forces, i.e., forces parallel to directions connecting 



2 Note that the vector u appearing in 13.48) can be replaced by any v £ O-u 



25 



the particle to its neighbors^ We will see later in Section [331 that the central-force 
decomposition is the only physically-meaningful partitioning of the force. 

The remaining question is how different potential energy extensions affect the 
force decomposition in ( 13.501 ). We have already established in (13.48b and d3.49l > that 
the force /™* is independent of the particular extension used. However, we show 
below that the individual terms in the decomposition, f a p, are not unique. These 
terms depend on the manner in which the potential energy, defined on the shape 
space, is extended to its neighborhood in the higher-dimensional Euclidean space. 

In order to construct different extensions, we use the geometric constraints that 
the distances have to satisfy in order for them to be embeddable in M? The nature 
of these constraints is studied in the field of distance geometry, which describes the 
geometry of sets of points in terms of the distances between them (see AppendixlBl). 
One of the main results of this theory, is that the constraints are given by Cayley- 
Menger determinants, which are related to the volume of a simplex formed by N 
points in an N — 1 dimensional space. 

For simplicity let us restrict our discussion to one dimension. It is easy to see that 
in one dimension the number of independent coordinates are N — 1 and for N > 2 
the number of interatomic distances exceeds the number of independent coordinates. 
Therefore, let the material system Ai consist of three point masses interacting in one 
dimension. The standard pair potential representation for this system, which is an 
extension of the potential energy to the higher-dimensional Euclidean space, is given 
by 

Vi„ t (Ci2, Ci3, C23) = Vi 2 (Ci 2 ) + V 13 (Ci 3 ) + V 23 (C2 3 ). (3.52) 
Since the calculation gets unwieldy, let us consider the special case when the particles 
are arranged to satisfy x\ < a; 2 < 23, for which 7*13 = ?*i 2 + r 23 . Using ( 13.50b . the 
internal force, f{ nt , evaluated at this configuration, is decomposed as 

dVi2 dV 13 

V 12 (ri 2 )+V 13 (r 13 ) 
/12 + /13. (3.53) 

23 The result that the force on a particle, modeled using any interatomic potential with a continuously 
differentiable extension, can be decomposed as sum of central forces may seem strange to some readers. 
This may be due to the common confusion in the literature of using the term "central-force models" to 
refer to simple pair potentials. In fact, we see that due to the invariance requirement stated on Page l21l 
all interatomic potentials (including those with explicit bond angle dependence) that can be expressed as 
a continuously differentiable function as described in the text, are central-force models. By this we mean 
that the force on any particle (say a) can be decomposed as a sum of terms, f a n, aligned with the vectors 
joining particle a with its neighbors and satisfying action and reaction. 

The difference between the general case and that of a pair potential is that for a pair potential, 
depends only on the distance r a p between the particles, whereas for a general potential, the dependence 
is on a larger set of distances, = ag'"^ ( ri2 > ri 3> • • ■ ' r (iV— l),JV*)i *- e -. II/q.sII depends on the 

environment of the "bond" between a and /3. For this reason, f a p for a pair potential is a property of 
particles a and /3 alone and can be physically interpreted as the "force exerted on particle a by particle 
0". Whereas, in the more general case of arbitrary interatomic potentials, the physical significance of the 
interatomic force is less clear and at best we can say that f a g is the "contribution to the force on particle 
a due to the presence of particle /3". 

24 We thank Ryan Elliott for suggesting this line of thinking. 



/i nt (n 2 ,ri3,r 23 ) = - 



"Vint 

dx\ 
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We now provide an alternate extension to the standard pair potential representation 
given in ( 13.521 ). The Cayley-Menger determinant corresponding to a cluster of three 
points (see (IB. 6b ) is identically equal to zero at every point on the shape space. This 
is because the shape space corresponds to a configuration of three collinear points, 
and the area of the triangle formed by three collinear points is zero. Thus, we have 

X(rL2, n.3, ^23) = (712 - r i3 - r 23 )(r 23 - r 12 ~ ^13) 
x (ri3 - r 2 3 - ?'i 2 )(?'i2 + ri 3 + r 23 ) 
= 0. (3.54) 

Using the identity in ( 13.54b . an alternate extension Vj£ t is constructed: 

V£t(Cl2, Cl3, C23) = V int (Cl2, Cl3, C23) + X(Cl2, Cl3, C23). 0.55) 

Note that Vj£ t is indeed an extension because from ( 13.54b it is clear that Vj£ t is equal 
to Vint at every point on the shape space of the system and it is continuously differen- 
tiable because x((i2, Ci3i C23). being a polynomial, is infinitely differentiable. Let us 
now see how the internal force, /} nt , for the special configuration considered in this 
example, is decomposed using the new extension: 

dVint _ dx 
dx\ dxi 

(/12 - 8ri 2 r 23 (ri 2 + r 23 )) + (/ 13 + 8ri 2 r 23 (ri 2 + r 23 )) 
A2+/13, (3.56) 

It is clear from (13.53b and ( 13.561 ) that the central-force decomposition is not the same 
for the two representations, i.e., /12 7^ /12 and /x3 7^ /13, however the force on 
particle 1, /{ nt , is the same in both cases as expected. 

It is very interesting to note that V{^ t is not a pair potential (based on the definition 
of a pair potential), but it is equivalent to a pair potential, i.e., it agrees with a pair 
potential on the shape space. Thus, the set of continuously differentiable extensions 
of a given interatomic potential function form an equivalence class. It is not clear at 
this stage if these equivalence classes can be fully expressed in terms of the Cayley- 
Menger determinant constraints. 

Although the above example is quite elementary, this process can be extended 
to any arbitrary number of particles in three dimensions. Any given potential can be 
altered to an equivalent potential by adding a function of the Cayley-Menger determi- 
nants corresponding to any cluster of 5 or 6 particles (see AppendixlBl. This function 
must be continuously differentiable and equal to zero when all of its arguments are 
zero. For example, a new representation in three dimensions can be constructed by 
adding a linear combination of the Cayley-Menger determinants: 

m 

V* nt = Vint(Cl2, ■ • ■ , C(JV-1)JV ) + ]TA fe Xfc, (3.57) 

k=l 



yint _ ""int 

' 1 dx\ 
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where there are m constraints defined by the Cayley-Menger determinants \k, and 
Afe are constants 1^1 

From this point on, we abuse our notation slightly, and write for any s = (ri2 , . . . , 
r(N-i)N) e S: 

dV m t dV m t , \ 

or a p dQap 

Also, we assume that there exists a continuously differentiable extension whenever 
we write f a p and sometimes refer to a continuously differentiable extension as an 
extension. 



Derivation of the pointwise stress tensor 



We now return to the differential equation in (13.36b for the potential part of the point- 
wise stress tensor. Substituting the force decomposition given in ( I3.501 l corresponding 
to a continuously differentiable extension, into ( 13.361 1, we obtain 



div x (T v (x, t) = ^2 (Wf a p \x a =x). (3.59) 

a,P 



On using the identity 



(fapW | x a = x) = l (fapW | x a = x, x f3 = y) dy, (3.60) 

Jm 3 



equation ( 13.59b takes the form 




= x,x g = y) dy. (3.61) 



We now note that, being anti-symmetric, the integrand in the right-hand side of the 
above equation satisfies all the necessary conditions for the application of Lemma 
IC. II given in Appendix A. Conditions (1) and (2) in Appendix A are satisfied through 
the regularity conditions on W. Therefore, using Lemma ICTl which was proved by 



25 Note that 13.57) has the same form as a Lagrangian with the A terms playing the role of Lagrange 
multipliers. For a static minimization problem, we seek to minimize V ; * lt , without violating the physical 
constraints relating the distances to each other. (This is equivalent to minimizing Vint with respect to the 
positions of particles.) Thus, the original constrained minimization of Vj n t is replaced by the problem of 
finding the saddle points of V;* t . 
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a 



8)Z 



Fig. 3.1 A schematic diagram helping to explain the vectors appealing in the pointwise potential stress 
expression in 43.62t . The bond a— ft is defined by the vector z. When s = 0, atom a is located at point x, 
and when s = 1, atom f5 is located at x. 



Noll in gU, we have 

<r v (x,t) (3.62) 
-fafiW \ x a = x + sz, xp = x — (1 — ds zdz 



2 7t ^ R3 J s =°\ dra P 



/ ( I !!!„ = a; + sz, xr = x — (1 — s)z) dsdz, 

Js=o \ dr afi I 



(3.63) 



where in passing to the second line, we have used ( 13.511 ) and the identity x a — xp = 
x + sz — [x — (1 — s)z] = z. For the special case of a pair potential, dVi n t/dr a p = 
V a p(r a p), and ( 13.621 ) reduces to the expression originally given in l48l . 

The expression for the potential part of the pointwise stress tensor in ( 13.62b is 
a general result applicable to all interatomic potentials. We make some important 
observations regarding this expressions below: 

1. Although the expression for er v appears complex, it is actually conceptually quite 
simple. er v at a point x is the superposition of the expectation values of the forces 
in all possible bonds passing through x. The variable z selects a bond length 
and direction and the variable s slides the bond through x from end to end (see 
Fig. ED. 

2. er v is symmetric. This is clear because the term z ® z is symmetric. Since the 
kinetic part of the stress in ( 13.35b is also symmetric, the conclusion is that the 
pointwise stress tensor is symmetric for all interatomic potentials. 

3. Since <r v depends on the nature of the force decomposition and different exten- 
sions of a given potential energy can result in different force decompositions, we 
conclude that the pointwise stress tensor is non-unique for all interatomic poten- 
tials (including the pair potential). We show in Sec tion [B31 that the difference due 
to any two pointwise stress tensors, resulting from different extensions for the in- 
teratomic potential energy, tends to zero as the volume of the domain over which 
these pointwise quantities are spatially averaged tends to infinity. Therefore, as 
expected, the macroscopic stress tensor, which is defined in the thermodynamic 



29 



limit (see footnote Q] on page[5]l, is always unique and is independent of the po- 
tential energy extension. 
4. Another source of non-uniqueness is that any expression of the form, er v + er, 
where div^ <r = 0, also satisfies the balance of linear momentum and is there- 
fore also a solution. We address this issue in Section I5T1 where we show that 
in the thermodynamic limit under equilibrium conditions, the spatially averaged 
counterpart to cr v converges to the virial stress derived in Section|2] 

The above results hinge on the use of the central-force decomposition in ( 13.501 ). 
One may wonder whether other non-central decompositions exist, and if yes, why 
are these discarded. This is discussed in the next section. 



3.5 Non-central-force decompositions and the strong law of action and reaction 

In the previous section, we showed that as a consequence of the invariance of the po- 
tential energy with respect to the Euclidean group, for any interatomic potential with 
a continuously differentiable extension, the force on a particle can always be repre- 
sented as a sum of central forces. In this section, we show that other non-central- 
force decompositions are possible, however that these violate the strong law of ac- 
tion and reaction, which we prove below, and therefore do not constitute physically- 
meaningful force decompositions. 

A proposal for a non-symmetric stress tensor for a three-body potential 

As an example, let us now consider the case of a three-body potential. For simplicity, 
we assume that the potential only has three-body terms and all particles are identical. 
Under these conditions, the internal potential energy is 

Vmt = ^2 V(x a ,xp,x^), (3.64) 

Q</3<7 

where V(x a , xp, x~) is the potential energy of an isolated cluster, {a,/3, 7}, and 
X) a,/3,7 represents a triple sum. We know that a central-force decomposition can 

a</3<7 

be obtained by following the procedure outlined in the previous section and that this 
leads to a symmetric pointwise stress tensor in ( 13.62b . Alternatively, a non-symmetric 
three-body stress tensor is derived as follows. To keep things simple, we derive the 
stress tensor for a system containing only three particles. Rewriting ( 13.64b for this 
case, we have 

3 

Vint = V(x 1 ,x 2 ,x 3 ) = ^2 4>a, (3-65) 

a=l 



where 



0a = ^{X1,X 2 ,X 3 ) 



(3.66) 
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is the potential energy assigned to particle a, equal to one-third of the total potential 
energy. Substituting ( 13.651 ) into ( 13.36b . we obtain 

div^ cr v (x,t) = ~y] {WV Xa 4>f} | x a = x) 

a,0 



^2 ( W ^<Ba^P I X a = X) - ^2 (W^Xc'fra \x a =x) 



ot±P 



(3.67) 



Since the cluster of three particles is isolated, the net force on the cluster due to 
internal interactions is zero. Therefore, from ( I3.651 l, we have 



= - ^2 v ^"- 



Using this relation, equation ( 13.67b simplifies to 

div x (T v (x,t) = - ^ (W(V Xa (/)p - V Xfl (f> a ) \x a = x) 



a, P 



Let 



/a/3 '■= ^Xf}4>a — V Xa <frp- 



(3.68) 



(3.69) 



(3.70) 



Now, using the identity 

(fapW | x a = x) = I (f a pW | x a = x, xp = y) dy , 



and the definition given in ( 13.701 ), equation d3.69t takes the form 

div^ er v (x, t) = V* / (Wf a p \ x a = x, xp = y) dy. 
T7> Jr 3 



(3.71) 



a7 tp 



Let us now study the definition of f a p given in ( 13.701 ). From ( 13.661 ) we have 



-V. 



dV dV 



dx a dxt 



The above equation suggests how the force /™* is decomposed. For example, f[ nt is 
decomposed as 



/i nt = /i2 + /is = ^(/i nt - fr) + x(/i nt - fT) 



Rearranging this relation gives 

/i nt + S'T + S'T = o, 

which is true since the cluster {1, 2, 3} is isolated. 



(3.73) 



(3.74) 
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1 1 




(a) (b) 



Fig. 3.2 (a) shows the force on each particle in a system consisting of 3 particles which interact through 
a 3-body potential given in 13.65) . Since the potential is derived from an energy decomposition, we have 
/} nt + f™* + f^ nt = 0. (b) shows the force decomposition of each /^ nt such that f a p = —fp a , but 
not necessarily parallel to the line joining particles a and /3, 



From (13.72b it is clear that / Q( g is anti-symmetric with respect to its arguments. 
Therefore, the integrand on the right-hand side of (13. 71b satisfies all the necessary 
conditions for the application of Lemma ICTI given in Appendix ICl Conditions (1) 
and (2) in Appendix[C]are satisfied through the regularity conditions on W. Therefore, 
using Lemma lCTl we have 

& v (x,t) = i V f f 

a 13 jR3 Js=0 

{-(Va^a - ^x a 4>p)W | x a = x + sz, x f} = x - (1 - s)z) ds ® zdz. 

(3.75) 

The stress <x v is non-symmetric in general because f a p, defined in ( 13.721 ), need not 
be parallel to the line joining particles a and /? as shown in Fig. 13.21 We therefore 
have two expressions for the stress for the same three-body potential. The symmetric 
expression in (13.62b and the non-symmetric expression in ( 13.75b . We show next that 
the non-central-force decomposition that led to the non-symmetric stress tensor is not 
physically meaningful since it violates the strong law of action and reaction. 

Weak and strong laws of action and reaction^ 

The following derivation hinges on the fact that in a material system the balance laws 
of linear and angular momentum must be satisfied for any part of the body. 

Consider a system of N particles with masses m a (a = 1, . . . , N). The total 
force on particle a is 

fa = + E f«f>> ( 176 > 



This derivation is due to Roger Fosdick 1201 . 
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where /° xt is the external force on particle a, and, as above, f a p is the contribution 
to the force on particle a due to the presence of particle (3. No assumptions are made 
regarding the terms f a p or the interatomic potential from which they are derived. 

A "part" p t of the system consists of K < N particles. We suppose Xq is a fixed 
point in space. Let F cxt (p t ) denote the total force on the part p t external to the part. 
Let M cxt (p t ; Xo) denote the total external moment on p t about Xq. Let L(p t ) be 
the linear momentum of the part p t and H(p t ; Xo) be the angular momentum of pt 
about Xq. 

We adopt the following balance laws, valid for all parts of the svstemF^I 

F cxt (p t ) = ^-(p t ), 0.77) 

J TT 

M^(p t ;x ) = —(p t ;x ). (3.78) 
at 

We now show that by applying these balance laws to particular parts of the system, 
that the strong law of action and reaction can be established. As a first observation, 
let pt consist of the single particle a. The external force and linear momentum for 

p t = {a} is 

f CXt (W) = /fW + E^«, (3-79) 

7 
7#« 

L({a}) = m a x a (t) (no sum). (3.80) 
The balance of linear momentum in ( 13.77b requires 

/a Xt + Yl f<*1 = m «*«' (3 ' 81) 

7 

The external moment and angular momentum of p t is 

M oxt ({a};a;o) = (x a (t) - x ) x (/^ xt + ^ f ai ) = m a (x a (t) - x ) x x a (t), 

(3.82) 



7 
7#a 



where we have used ( 13. 8U . and 

H({a}; xo) = (x a — x ) x m a x a (t). (3.83) 
The balance of angular momentum in (13.78b is satisfied identically, since 

m a (x a (t) - x ) x x a (t) = ^ [(x a (t) - x ) x m a x a (t)} 

= x a (t) x m a x a (t) + (x a (t) — Xo) x m a x a (t) 
= m a (x a (t) - x ) x x a (t). 



27 The view that the balance of linear momentum and the balance of angular momentum are funda- 
mental laws of mechanics lies at the basis of continuum mechanics. See, for example, Truesdell's article 
"Whence the Law of Moment and Momentum?" in 1611 . 
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As a second observation, let p t consist of the union of the two particles a and (3. 
The external force and linear momentum are 

F cxt ({a,/3}) = /r + E (/»7 + //»t)» ( 3 - 84 > 

L({a,/3}) = m a x a + mpxp. (3.85) 
The balance of linear momentum in ( 13.77b requires 

/ c Xt + / cxt + ^ (fa + fa) = m a x a + mpxp. (3.86) 

Subtracting (13.8 lb for particles a and (3 gives 

E (/or + ffh) ~ E /«t - E #7 = 0= (3-87) 

7 7 7 

from which 

/a/5 + //3a = 0. (3.88) 

This relation is called the weak law of action and reaction f22\ . It shows that f a p = 
—fpa, but does not guarantee that f a p lies along the line connecting particles a and 

Next, the external moment and angular momentum of p t is 

= (x a -x )x(f^ + f^) + (xp-x )x(f^ + Y, fa) 

7 7 

= (x a - x ) x (m a x a ~ fa) + (xp - x ) x (mpxp ~ fa), (3.89) 
where we have used ( 13.811 ), and 

H({a, P}; x ) = (x a - x ) x m a x a + (xp - x ) x mpxp. (3.90) 
The balance of angular momentum in ( 13.78b requires 

(x a - x ) x (m a x a - fa) + (xp - x ) x {mpxp - fa) 
= [(x a - x ) x m a x a + (xp - x ) x mpxp] 

= x a x m a x a + (x a — xo) x m a x a + xp x mpxp + (xp — Xq) x mpxp 

= (x a - xo) x m a x a + (xp - xq) x mpxp, (3.91) 

which simplifies to 

(X a - X ) X fa + (xp - Xq) X fa = 0, 



34 



and, after using ( 13.88b . we obtain 

(x a - xp) x f aP = 0. (3.92) 

This shows that f a p must be parallel to the line joining particles a and /3. This is 
the strong law of action and reaction. We have shown that this law must hold for 
any force decomposition, in order for the balance of linear and angular momentum to 
hold for any subset of a system of particles. 

The possibility of non- symmetric stress 

Based on the proof given above for the strong law of action and reaction, we argue that 
only force decompositions that satisfy the strong law of action and reaction provide 
a physically-meaningful definition for f a p. For example, the definition in ( 13.701 ) is 
not physical because if it were used to compute the external moment acting on a sub- 
system of particles, as is done above, the balance of angular momentum would be 
violated. For this reason, this decomposition and the corresponding non-symmetric 
stress in (13.75b are discarded. The conclusion is that a pointwise stress tensor for a 
discrete system of point masses without internal structure has to be symmetric. 

In the next section, we discuss the possibility of expanding the class of solutions 
resulting from Irving-Kirkwood-Noll procedure in a way that makes it possible to 
obtain non-symmetric stress tensors for systems where the point particles have inter- 
nal structure. This involves a relaxation of the assumption that the "bonds" connecting 
two particles are necessarily straight. 

3.6 Generalized non-symmetric pointwise stress for particles with internal structure 

In Section 13.41 we saw that the Irving-Kirkwood-Noll procedure, when applied to 
multi-body potentials, results in a symmetric non-unique pointwise stress tensor. We 
now seek to find additional solutions to d3.611 l. which are not obtained using the 
standard Irving-Kirkwood-Noll procedure. In arriving at ( 13.62b using Lemma |CT1 
we can see that the contribution to the potential part of the stress at position x is 
due to all possible bonds, assumed to be straight lines, that pass through x. The 
question that naturally arises is to what extent can this assumption be weakened. In 
other words, can Lemma lCTl be generalized in a suitable manner so that non-straight 
bonds can be accommodated? Such a possibility was first discussed by Schofield and 
Henderson 1521 . who used the Irving and Kirkwood approach with a series expansion 
of the Dirac-delta distribution. It will be shown in this section, using Noll's more 
rigorous approach, that solutions giving rise to non-straight bonds are possible. 

From a physical standpoint, non-straight bonds are possible for systems with in- 
ternal degrees of freedom. An example would be the dipole-dipole interactions be- 
tween water molecules resulting from the electrical dipole of each molecule. The 
possibility of internal degrees of freedom was already raised by Kirkwood in his 
1946 paper |29l . The idea is to relate the shape of the non-straight bonds to the ad- 
ditional physics associated with the internal degrees of freedom. This issue will be 
further explored in future work. For now, we only investigate the possible existence 
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of additional solutions other than that given by (13.62b . We begin by describing the 
shape of a bond in a more precise way through the following definition. 

Definition 2 The "path of interaction" between any two interacting particles a and 
(5 is the unique contour that connects a and f3, such that there is a non-zero contri- 
bution to the potential part of the pointwise stress, <x v , at any point on this contour. 

In this section, the terms bond and path of interaction are used synonymously. 
Therefore, for the case of the pointwise stress tensor in ( 13.621 ). this path of interac- 
tion is given by the straight line joining a and (5. It is shown in Appendix ICl that 
under certain restrictions on the path of interaction, Lemma ICTI can be generalized 
to Lemma |Q21 given in Appendix ICl Roughly speaking these restrictions are given 
by the following conditions: 

1 . The shape of the bond connecting particles a and j3 only depends on the distance 
between a and j3. o 

2. For any two pairs of particles (a, j3) and (7, 6) separated by the same distance, 
the bonds a — (3 and 7 — S are related by a rigid body motion. In addition, if 

then this rigid body motion involves only translation. 

From conditionQ] it is clear that the shape of the bonds can be described by contours 
Yi : [0, 1] -> R 3 , for I > 0, with Tj(0) = (0, 0, 0), T ; (l) = (I, 0, 0) along with some 
mild restrictions. Hence, defining the contours T; for / > is equivalent to defining 
the paths of interaction between any two points in R 3 . For a precise definition of Ti 
and the paths of interaction see Appendix [C] Since all the necessary conditions for 
the application of Lemma IC2l are same as those for Lemma ICTI we can use Lemma 
IC.2l in the Irving-Kirkwood-Noll procedure instead of Lemma ICTI In doing so, we 
arrive at a definition for the generalized pointwise stress tensor cr^, for given paths 
of interaction with the above mentioned properties. This is given by 

a°(x,t) = 

\ 5^ / / (f°>pW \x a = x± + sz, xp =x±- (1- s)z) <g) Q*TL|| (s) ds dz, 

a p ,/r3 Js=0 

(3.93) 

where f a p is defined in ( 13.511 ). and x±(s,x, z) = x — sz — Q^Tiuii (s), Q z E 
SO(3). Here, x± represents the projection of x onto the line joining the endpoints, 
x a = x± + sz and xp = x± — (1 — s)z, of the path of interaction being considered. 
Q z represents the rotation part of the rigid body motion described in condition|2]that 
maps the contour T|| a . Q _ a , a || to the path of interaction that connects x a and xp. 

Equation (13 ,93b is a general expression for the potential part of the pointwise 
stress tensor, of which ( 13.621 ) is a special case. We discuss several key features of this 
expression below: 

28 This is in essence a constitutive postulate similar to the assumption in pair potentials that the energy 
depends on only the distance between particles. A more general dependence of the shape of the path of 
interaction on the environment of a pair of particles might be possible, but is not pursued here. See also 
Definition |3]in Appendix Icl and following discussion. 
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1 . Equation (13.93b is unique only up to given paths of interaction for a given poten- 
tial energy extension. It is a more general result than (13.621 i. since ( 13.621 ) can be 
obtained from (13.93b by assuming that the path of interaction between any two 
points is the straight line connecting them. For this special case it is easy to see 
that x± = x and Q z Y^ z ^ (s) = —z. 

2. cr^ is in general non-symmetric, whereas the stress obtained through the standard 
Irving-Kirkwood-Noll procedure is always symmetric for any multi-body poten- 
tial with an extension. Since the kinetic part of the stress tensor er^ (see ( 13.35b ) 
is symmetric, it follows that the total pointwise stress tensor obtained from the 
generalized stress tensor is usually non-symmetric. Therefore, under the present 
setting, the balance of angular momentum is satisfied only through the presence 
of couple stresses. This suggests that non-straight bonds might correspond to sys- 
tems with particles having internal degrees of freedom. 

3. Since both ( 13.62b and ( 13.931 ) are valid definitions for the potential part of the 
pointwise stress, the question of which one to choose depends on the presence of 
internal degrees of freedom in each particle. In the absence of internal degrees 
of freedom only straight bonds are possible due to symmetry. The issue of the 
pointwise stress being unique only up to a divergence-free tensor-valued function 
is partially addressed here, since the expression given by the difference between 
the two definitions is divergence-free. 

4. The expression in ( 13.93b is very similar to ( 13.621 ). The pointwise stress at x is 
a superposition of the expectations of the force of all possible bonds/paths of 
interaction passing through x. The vector z selects an orientation and the size of 
the vector connecting the two ends of the bond, and s slides it through x from 
end to end as shown in Fig. 13.31 



3.7 Definition of the pointwise traction vector 

In this section, we derive the formula for the pointwise traction vector t(x,n;t) 
defined on the surface passing through x, with normal n at time t. The following 
derivation is based on [48 1 and it can be easily extended to curved paths of interaction. 
As usual, let Ai denote our material system. Let fl C R 3 be a domain in three- 
dimensional space with continuously differentiable surface S, representing a part of 
the body. By this definition, each of the N point masses described by Ai either belong 
to J? or in the space surrounding J?, denoted by $7 C . Let / denote the force exerted 
by the particles in J? c on particles in Q. We note that in continuum mechanics / is 
related to t by 



where n(x) is the outer normal at x 6 S. Using the Cauchy relation, t(x, n, t) = 
cr(x, t)n, we obtain 




(3.94) 




(3.95) 
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(a) 

3 




(C) 

Fig. 3.3 A schematic diagram helping to explain the vectors appearing in the inner integral of 0.93) for a 
given point x. The integral in 13.93) in an integral over all possible paths of interaction that pass through 
the point x. The inner integral with respect to s, with z fixed, is an integral over those paths, where z is 
the vector joining its endpoints. Frame (a) shows a path of interaction when s = 0. As s is increased the 
path "slides" through x. Frame (b) shows the path for an arbitrary s in the interval [0, 1]. The end points 
are represented by + sz and x± — (1 — s)z. Frame (c) shows the position of the path for s = 1. 



Now, note that the net force exerted by i? c on f2 due to particle interaction, denoted 
by f v (t), is given by 

/v(*) = X] / / (fapW I x a = u,x p = v) dudv, (3.96) 



3K 



where f a p is defined in ( I3.51l l. Since the integrand in ( I3.96l l satisfies all the conditions 
for the application of the lemmas in Appendix [C] we can now use a special case of 
Lemma lC31 bv restricting to straight bonds0 We therefore have, 



Mt) 



2 



S JR 3 Js=0 



/ ( — fapW | x a = x + sz, xp = x — (1 — s)z) (z ■ n) ds dz dS(x). 



(3.97) 



We now note that f v in ( 13.971 ) exactly satisfies 

Mt) = [ <r v ndS{x), (3.98) 
Js 

where <r v is given by ( 13.621 ). It is therefore clear that f v describes the potential part of 
the interaction force /. Hence, it is natural to assign a potential part of the pointwise 
traction vector, t v to / v , given by 

t v (x, n; t) := cr v n 

= \ / / (-fapW | x a = x + sz, Xp = x - (1 - s)z) (z ■ n) ds dz. 

(3.99) 

The above formula is conceptually quite simple. It gives the measure of the force 
per unit area of all the bonds that cross the surface, where the force is calculated 
with respect to a surface measure (see ( 13.97b ). Using this viewpoint, we motivate 
the definitions for the macroscopic traction vector and the stress tensor, when we 
incorporate spatial averaging in the next section. 

It is now natural to assign the kinetic contribution to the force across the surface 
to the kinetic part of the pointwise stress tensor. Subtracting ( 13.981 ) from ( 13.951 ), we 
obtain the kinetic contribution to the force across a surface, 

fk(t) := / (<r — <r v )ndS(x) = / er k n dS(x). 
Js Js 

Therefore the kinetic contribution to the pointwise traction vector is given by 
tk(x, n; t) := er k n 

= -J2 ™ a «'« cl • n)W \x a =x). (3.100) 

a 

Finally, we note that the definitions of t v and ty are functions of x and n alone. 
Hence, this result is related to the work of Fosdick and Virga I2TI . who give a varia- 
tional proof for the stress theorem of Cauchy in the continuum version. In that work 



'Specifically, for straight bonds, we set: x±_ = x and Q z r! ,, (s) 
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the traction vector is allowed to depend on the unit normal and the surface gradient 
and is shown to be independent of the surface gradient. 

The fields defined and derived in this section are pointwise quantities. In the next 
section, expressions for macroscopic fields are obtained by spatially averaging the 
pointwise fields over an appropriate macroscopic domain. 



4 Spatial averaging 

In the previous section, the Irving-Kirkwood-Noll procedure was used to construct 
pointwise fields from the underlying discrete microscopic system using the principles 
of classical statistical mechanics. Although the resulting fields resemble the contin- 
uum mechanics fields and satisfy the continuum conservation equations, they are 
not macroscopic continuum fields. For example, the pointwise stress field in ( 13.62b . 
at sufficiently low temperature, will be highly non-uniform, exhibiting a criss-cross 
pattern with higher stresses along bond directions, even when macroscopically the 
material is nominally under uniform or even zero stress. 

To measure the fields derived in the previous section in an experiment, one needs 
a probe which can extract data only from a single point of interest in space. Since 
this is not possible practically, there is no way we can correlate the experimental data 
with theoretical predictions. Therefore a true macroscopic quantity is by necessity 
an average over some spatial region surrounding the continuum point where it is 
nominally definedQ Thus, if f{x,t) is an Irving-Kirkwood-Noll pointwise field, 
such as density or stress, the corresponding macroscopic field f w (x, t) is given by 



where w(r) is a weighting function representing the properties of the probe and its 
lengthscale. 

The important thing to note is that due to the linearity of the phase averaging in 
the Irving-Kirkwood-Noll procedure, the averaged macroscopic function f w (x,t) 
satisfies the same balance equations as does the pointwise measure f(x, t). 

Weighting function 

The weighting function w(r) is an R + -valued function with compact support so that 
w(r) = for ||r|| > A, where A is a microstructural lengthscale. The weighting 
function has units of volume -1 and must satisfy the normalization condition 



We do not include time averaging, because this is indirectly performed due to the presence of W. 
The reasoning for this comes from the frequentist's interpretation of probability, wherein the probability 
of a state is equal to the fraction of the total time spent by the system in that state. 




(4.1) 




(4.2) 
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0.5 1 1.5 

r/r w 

Fig. 4.1 Three weighting functions for spatial averaging: uniform weighting (solid line) in j4.3t : Gaussian 
weighting (dashed line) in i4A\ ; Quartic spline weighting (dash-dot line) in i4.5\ . Note that the areas under 
the curves are not equal because the normalization in )4.2t is according to volume. 



This condition ensures that the correct macroscopic stress is obtained when the point- 
wise stress is uniform. For a spherically-symmetric distribution, w(r) = w(r), where 
r = \\r\\. The normalization condition in this case is 

w(r)Airr 2 dr = 1. 

The simplest choice for w(r) is a spherically-symmetric uniform distribution over a 
specified radius r w , given by 



> (r) = /l/K.ifr< 



i l/v w ii r s r w , 
\ otherwise, 



where V w = jirr^ is the volume of the sphere. This function is discontinuous at 
r = r w . If this is a concern, a "mollifying" function that smoothly takes w(r) to zero 
at r w over some desired range can be added H31 F*l Another possible choice for w(r) 
is a Gaussian function I 



w(r) = 7T ^ r J cxp [- r yri\ . (4.4) 

This function does not have compact support. However it decays rapidly with distance 
so that a numerical cutoff can be imposed where its value drops below a specified 
tolerance. Another possibility is a quartic spline used in meshless method applications 
(where it is called a kernel function Q), 

w(r) = I 167 "» r » A ~ ' (4.5) 

l_ otherwise. 

This spline has the advantage that it goes smoothly to zero at r = r w , i.e., w(r w ) = 
0, w'(r w ) = 0, and w"(r w ) = 0. Fig. 14. II shows the plots of the three weighting 
functions given above. 



31 An example of a mollifying function is given later in equation j6.13t . 
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4. 1 Spatial averaging and macroscopic fields 



Continuum fields such as density and momentum density fields are defined using 
( 14. Il l as the ensemble average via the probability density function W, followed by a 
spatial average via the weight function w as follows: 



p w (x, t) := y m a / w(y - x)(W | x a = y) dy, 
p w (x,t) := Y]m Q / w(y - x)(Wv a \ x a = y) dy. 



(4.6) 
(4.7) 



It is straightforward to show that using definitions, (14.61 l and (14.7b . the macroscopic 
version of the generalized pointwise stress tensor given by ( 13.931 l divides into poten- 
tial and kinetic parts as, 



Om) = ^E [ w (v~ x ) [ I 

1 „ o JS? J«3 Js=0 



(fapW | x a = y± + sz, xp = y± - (1 - s)z) ® QzT^ (s) ds dz dy, 

(4.8) 

where f a p is defined in (13.51b , y± = y — sz — QzY\\ z n (s), Q z G SO(3), and 

cr Wik (x,t) = - J2 I w(y- x)m a ((v™ 1 ® < el ) W | * a = 1/} dy . (4.9) 

We now intend to express the potential part of stress in a more convenient form. This 
is done by two consecutive changes of variables. Under the assumption that Q z and 
f\\z\\ are differentiable with respect to z and ||z||, respectively, the Jacobian of the 
transformation (s, y, z) M> (s, y±, z) is unity. Therefore, 



v(rc,t) = \y\ ( w{y-x) ( [ 



a,/9 

| a; Q = y ± + sz,x^ = y± - (1 - s)z) ® Qz^\\ z \\ dsdzdy±, 

(4.10) 

where y = y(s, y±, z). A second change of variables is introduced as follows 

y±+,sz = u, y±-(l-s)z = v, (4.11) 

which implies, 



z = u — v, y± = (1 — s)u + sv. 
The Jacobian of the transformation is 



J = det 



V u z W v z 



= det 



J —I 

(1 - s)I si 



= 1. 



(4.12) 



(4.13) 
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w((l — s)u + sv — x) 




Fig. 4.2 The bond function b(x: u, v) is the integral of the weighting function centered at x along the 
line connecting points u and v. The graph shows the result for a quartic spline weighting function. The 
bond function is the area under the curve. 



where 



Using ( 14.11b , ( 14.12b and ( 14.13b to rewrite (14.10b . we obtain 

<r Wj y(x,t) = - ^ / (-fapW | x a = u,xp = v) <g) b{x;u,v)dudv, 

a, 8 J R3xR3 

(4.14) 

b{x;u,v) := - / u>(y - a;)Q u _^r,L_^|| (4.15) 
is called the bond vector, with 

y(s,u,v) = y(s,yx(s,u,v),z(u,v)). 
For the special case of straight bonds, we have 

y = (1 - s)u + sv and Q U - V T^ U _ V ^ (s) = -(« - v). 
Therefore the bond vector simplifies to 

b(x; m, v) = (u — v) / — s)u + sv — x) ds 

Js=0 

= (u — v)b(x;u,v), (4.16) 

where b(x; u, v) is commonly referred to as the bond function. The geometrical sig- 
nificance of the bond function is explained in Fig. 14. 21 

For the special case of straight bonds, equation ( 14.141 ) simplifies to 

tr w , v (x,t) = -J2 / (-fapW I x a = u,xp = v)®(u-v)b(x;u,v) dudv. 

(4.17) 

The expressions for the potential and kinetic parts of the spatially-averaged stress 
tensor in equations d4.9b and ( 14.17b are our main result and constitute the general 
definitions for the macroscopic stress computed for a discrete system. It will be shown 
in Section[5]that these relations reduce to the Hardy stress tensor l23l under suitable 
approximations. The issue of the uniqueness of the stress tensor (in the sense that any 
divergence-free field can be added to it) is deferred to Section[ 
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4.2 Comparison with the Murdoch-Hardy procedure 

An alternative procedure of defining continuum fields to the one described above, due 
to Murdoch H43H46ll47l and Hardy l23l . only involves spatial averaging. We refer to 
this approach as the Murdoch-Hardy procedure. Under the Murdoch-Hardy proce- 
dure, continuum fields are defined as direct spatial averages of microscopic variables 
without incorporating statistical mechanics ideas. Therefore, the Murdoch-Hardy 
procedure is purely deterministic in nature. For example, the density and momen- 
tum density fields at a particular instant of time, corresponding to a given weighting 
function w, are defined as 



respectively, where x a and v a are deterministic quantities. We denote spatially- 
averaged variables obtained from the Murdoch-Hardy procedure with a superposed 
tilde to distinguish them from quantities obtained in Section |4~T1 Equation (14. 1 8b is 
used to "smear" a discrete system to form a continuum. The reasoning for abandon- 
ing statistical mechanics is the lack of knowledge of the ensemble of the system as 
explained by Murdoch and Bedeaux in 11461 : 

Physical interpretations of any given ensemble average clearly depends 
on the definition of the ensemble . . . for example, if a container is filled to a 
given level with water and then poured onto a surface, the lack of precision 
with which the pouring is effected may result in many different macroscopic 
flows. Here no single description is available within deterministic continuum 
mechanics: in this case the ensemble (defined in terms of the water molecules 
and limited knowledge of how the pouring takes place) relations involve aver- 
ages associated with all possible flows. Clearly relations involving ensemble 
averages are associated with a much greater variety of behavior than is de- 
scribable in terms of deterministic continuum mechanics. 

We share the same concern regarding the ambiguity in the definition of an ensemble. 
For example, in an experiment where an austenite-martensite phase transformation 
occurs, the resulting micro-structure consists of a complex spatial configuration of 
martensitic variants, and this depends largely on the microscopic details of the sys- 
tem, such as cracks, lattice defects, etc. Therefore, in this case, macroscopic variables 
cannot completely describe the ensemble of interest. To avoid this difficulty, Mur- 
doch proposes a time average in place of ensemble average. Nevertheless, it should 
be noted that from classical statistical mechanics, the ensemble of interest and its 
corresponding distribution exists in principle. Therefore the framework described in 
Section l4~T1 is a correct framework in which to phrase the problem. A practical calcu- 
lation can then be performed, for example, by replacing the ensemble averages with 
time averages in a molecular dynamics calculation (see Section|5]l- We stress the im- 
portance of writing a continuum field variable as an ensemble average followed by 
spatial average, rather than a spatial average followed by a time average, as is done 




(4.18a) 



a 




(4.18b) 



a 
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in the Murdoch-Hardy procedure, because it helps to give a unified picture of all the 
previous definitions for continuum fields and stress in particular. This is discussed in 
the next section. 

It is interesting to note that by relaxing the connection with statistical mechanics, 
the Murdoch-Hardy procedure allows for a much wider class of definitions for the 
stress tensor B31 in addition to the non-uniqueness characterized so far, due to the 
presence of multiple extensions for the potential energy and allowing non-straight 
bonds. In this section we intend to systematize this procedure. The source of non- 
uniqueness resulting from multiple definitions of the stress tensor is studied, thus 
helping us to identify a much larger class of possible definitions. In this new system- 
atic approach, the steps involved in the Murdoch-Hardy procedure are as follows: 

1. Develop a continuum system by smearing out the discrete system using ( 14. 181 ). 

2. Introduce a non-local constitutive law for the continuum that is consistent with 
the discrete version of force balance given later in ( 14. 191 ). 

3. For each constitutive law, define a stress tensor, which satisfies the equation of 
motion for the continuum. 

To understand the above three steps, we explore the Murdoch-Hardy procedure 
in more detail. The continuity equation is satisfied in a trivial way l44l . We now look 
at the equation of motion. 

Equation of motion 

The motion of particle a is governed by Newton's second law, 

^2f^(t) + b d a (t) = m a v a (t), (4.19) 

where := f a p(u{t)), f a p{u) are the terms in the central-force decomposition 

obtained from a multi-body potential with an extension (see Section l3~4l and b^(t) 
is defined as 

K(t) := -V^VextOciW, • ■ • , x N {t)). 

The superscript "d" in (14.19b and the above equation are used to the stress the fact that 
the quantities are deterministic in nature. Equation ( 14.19b is a force balance equation 
for the discrete system. We now design an analogous force balance equation for the 
smeared continuum defined by (14.181 ), such that ( 14.191 ) always holds. 

For the sake of simplicity in notation, from here onwards we use f a p to denote 
both f a p{u) and fap(t)> whenever it is clear from the context. The same goes with 
the usage of b a (t) for b^(t). 

Force balance for the smeared continuum 

Multiplying (14.19b by w(xi — x) and summing over all particles, we have 

fw + K, m Q v a w(x a - x), (4.20) 
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where 

f w (x,t) :=^2f a p(t)w(x a (t)-x), (4.21) 



ot,f) 



b w (x,t) :=^2bp(t)w(x p (t) - x). (4.22) 



To arrive at a form similar to the equation of motion of continuum mechanics given 
in ( 13.23b . equation (I4.20l i is rewritten as 

d 

fw + b w = — m a w(x a - x)v a - ^ m a v a CVw(x a - x) ■ v a ) 



dt 

dp w 
dt 



diva, m a w(x a — x)v a (8 v a . (4.23) 



Similar to ( 13.16b . we define the continuum velocity as 

v w {x,t):=¥^r, (4.24) 

and the relative velocity of a particle with respect to the continuum velocity as 

v™\x,t) := v a {t) - v w (x,t). (4.25) 
Using ( 14.181 1 and ( 14.25b . we obtain 

^m a v' c l fw(x a (t) -x)= p w (x,t) - p w (x,t)v w (x,t) 

OC 

= o, 

the last equality being true by the definition ( 14.24b . From (14.25b and the above equa- 
tion, it follows that 

rn a w(x a - x)v a <g> v a = ^ m a w(x a - x)v T ° l ® i>„ cl + ® 

a a 

Substituting this into ( 14.23b and rearranging, we have 

fw - div x ^ m a (v™ 1 (g> ■y™ 1 )w(a; Q - as) + b lu = + div x (p w v w ® 

a 

Comparing the above equation with the equation of motion, (13.23b . we have 

div x a w (x,t) = f w (x,t) - div x ^2m a (v^ 1 ® ^ el )w(x a - aj), (4.26) 
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where cr w is the stress tensor corresponding to the weighting function w. From ( 14.261 ) 
it is clear that the kinetic part and the potential part of the stress tensor, cr^k and & w>v , 
respectively, are given by 

v w jt(x,t) = -^m a (v I f®v I f)w(x a -x), (4.27a) 

a 

div x d- WtV (x,t) = f w (x,t). (4.27b) 

Any solution to ( I4.27bt is a valid candidate for the definition of a w>v . Murdoch B31I 
proposes several possible candidates, and highlights the possibility of having multiple 
definitions. To understand the connection between the different possible definitions, 
we look back at ( 14.20b and ( 14.211 ). Equation ( 14.201 ) is a force balance equation for 
any "continuum particle" at x, and f w , defined in ( 14.21b . is the force per unit volume 
acting on it. It is not immediately clear from ( 14.211 ) how two continuum particles at 
positions x and y interact with each other. This interaction can be given by a non- 
local constitutive law. The main idea is to recast ( 14.211 ) as 

f w {x,t)= [ g(x,y,t)dy, (4.28) 

for some g(x, y, t), which we call the generator of the non-local constitutive law. 
This function describes the interaction between the continuum particles at x and y . To 
satisfy Newton's third law, we also need g to be anti-symmetric with respect to its ar- 
guments x and y. Unfortunately the representation given in (14.28b is not unique and, 
since every choice of g leads to a different stress definition, this is one of the sources 
of non-uniqueness in the definition for the stress tensor in the Murdoch-Hardy pro- 
cedure. We describe two different constitutive laws, which lead to the Hardy stress 
and the doubly-averaged stress (DA stressj^ ED- 

For the case of Hardy stress, the generator g H is given by the equation 

g H (x, y, t) = f a (iw{x a - x)5(xp - x a + x - y), (4.29) 

a,0 

where 6 denotes the Dirac delta distribution. 

The generator g° for the DA stress is given by 

Q^ix, y,t) = ^2 faf}w{x a - x)w(xp - y). (4.30) 

Fig. l4.3l shows the interaction between two continuum particles, with positions x 
and y, that are in a neighborhood of two interacting particles a and f3 respectively 
and not in a neighborhood of any other particle in the system. In this setup, it is 
clear from the generator for the Hardy stress, given in (14.29b , that two continuum 
particles at x and y interact only when y — xp = x — x a , as shown in Fig. |4.3(a)| 



32 Murdoch 1451 refers to this stress as "Noll's choice". To avoid confusion with the stress derived 
through the Irving-Kirkwood-Noll procedure in Sectionf5] we name it the "DA stress". 
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(a) 



(b) 



9{x,y,t) 



(C) 



Fig. 4.3 A continuum particle x interacts with: (a) only that continuum particle at y, which is identically 
oriented to xp as x is oriented to x a , when the interaction is given by g H ; (b) any continuum particle 
in the shaded region, when the interaction is given by g° ; (c) any continuum particle on the shaded line, 
when the interaction is given by g HD . 



On the other hand, there is no such restriction on the generator for the DA stress, 
described by d430b (see Fig. |4.3(b)| >H| Although at this point there is no systematic 
way of suggesting additional possible generators, we can suggest a third generator, 
g HD , which has properties that lie in between g H and g° . As shown in Fig. 14.3(c)] 
when interaction is governed by </ HD , a continuum particle x interacts with y only 
when y lies on the line passing through x and parallel to x a — xp. In all the three 
cases, the interaction force is always directed along the vector x a — xp. Therefore 
by (14.28b we have three different integral representations for f w with generators g° , 
S H ,9 HD 



Now, in each of the integral representations of f w given by ( 14.291 ) and ( 14.301 ). the 
integrand satisfies all the necessary conditions for the application of Lemma ICTI or 
Lemma lC2l in Appendix lcf^l For instance, using Lemma lCTTI we obtain an expression 
for the potential part of the stress tensor, given by 



(x,t) = -- / g(x + sz,x-(l-s)z,t)ds 
z Jm 3 Us=o 



z dz. 



(4.31) 



33 Note that for the DA stress, the force between two continuum particles at x and y is not parallel to 
x — y in general. This is not a violation of the strong law of action and reaction, because the strong law 
only applies to discrete systems. It has been used in this derivation by requiring that f a = —fp a and 

fa/3 X (x a ~ Xp) = 0. 

34 The integrand should be continuously differentiable for Noll's lemma to be applicable. Although g H 
is not continuously differentiable due to the presence of the Dirac delta distribution, this does not hinder 
us from applying the lemma since we can replace the Dirac delta distribution by an appropriate infinitely 
differentiable delta sequence and take a limit. See Appendix iDlf or a rigorous derivation of this. 
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Substituting (I4.29l l into ( 14.3 It . we have the potential part of the Hardy stress: 

»i 

f a pw(x a — x — sz)5(x(s — x a + z) ds ® zdz 



w.v 



2 ^ 

It, f 



a, 



[~f a pw((l - s)x a + sxp - x) ® (x Q - xp)] ds. (4.32) 



Substituting (14. 30t into ( 14. 3U we have the potential part of the DA stress: 

- J2 I [ [-f a pw(x a -x-sz)w(xp-x + (l-s)z)(g>z]dsdz, 
z „ a Jzm 3 Js=0 



-DA 
w ,v 



a,0 

(4.33) 

which was derived by Murdoch flTl . The conclusion is that the non-uniqueness of 
the generator in the systematized Murdoch-Hardy procedure leads to a non-unique 
definition for the stress tensor. Further sources of non-uniqueness can by introduced 
by having a different force decomposition corresponding to a different potential ex- 
tension, or using curved paths of interaction instead of straight bonds and applying 
Lemma IC.2I which is a generalization of Lemma IC.ll in Appendix [C] We do not 
pursue this generalization further here. 

It is important to point out that the systematized Murdoch-Hardy procedure pre- 
sented here does not describe all possible solutions to (I4.27bt . An example of a so- 
lution that cannot be obtained via our systematized Murdoch-Hardy procedure is the 
following definition suggested in l45l : 



>v (x,t) := ^2 fafi ® i X 



t )a(\\x 



II), 



(4.34) 



a,/3 



where a(u) :— J " s 2 w(s) ds. As is pointed out in B31 . the expansion in ( 14.34b is 
not a physically-relevant definition for stress due to the following test case. Consider 
a stationary deformed body at zero temperature (i.e., where the particles occupy fixed 
positions without vibrating). In this case, the net force acting on any particle is zero. 
Since f a p is the only term in the summand of d4.341 i which depends on j3, (14.341 i is 
equivalent to 

oCvOM) : = ®( x ~ x a)a(\\x - Xo,||). (4.35) 

a Si 

In our case, f a p = f a = 0, for each particle a in the interior of the body 
which is considerably away from the surface compared to the interatomic distance. 
Hence, the only non-zero contribution to the stress is due to those particles close to 
the surface on which the net force due to other particles is non-zero. Moreover, al- 
though d(u) decays to zero as u increases, a(u) ^ for all u ^ 0. Thus, there is a 
non-zero stress at every point x G M 3 (even outside the body!) due to particles close 
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to the surface of the body, which is obviously not physically reasonable. Neverthe- 
less, er^ v is mathematically still a valid definition since it satisfies the force balance 
equation B31 . However, it cannot be derived using the systematized Murdoch-Hardy 
procedure proposed here. Thus, the systematized Murdoch-Hardy procedure does 
not lead to all possible definitions that satisfy the force balance equation. 

Finally, it is also worth noting that the fact that all balance laws are satisfied 
under the Murdoch-Hardy procedure should not come as a surprise, since w in ( 14.181 ) 
serves the same purpose as W in ( 13.13b . In this view, w is seen as a function defined 
on a phase space, although one that does not evolve according to a flow described 
by Hamilton's equations of motion, but still satisfies (13.10bP 5 l The corresponding 



flow in phase space is described as follows. Continuing with our notation introduced 
(ED, let S(0) = (a; s (0),« 8 (0)) = (xf(0), . . . ,afy(0),«!(0), . . . , v%(0)) denote 
any arbitrary point in phase space. We add a superscript "s" to stress the fact that an 
element in phase space is stochastic in nature. Consider the flow in phase space given 
by the mapping 

S(0) = (a; s (0), v s (0)) ^(x s (0) + x{t) - x{0),v s (0) + v{t) - v(0)) 

= (x s (t),v s (t)) = 3(t), (4.36) 



where the quantities x (t) = (xi(t),--- ,xpj(t)),v(t) = (vx(t),--- , v^(t)) denote 
the position and velocity of the particle and these are assumed to be known. (Typically 
these quantities are obtained from a molecular dynamics simulation.) Therefore the 
Murdoch-Hardy procedure can be interpreted as a probabilistic model constructed 
from the data, x(t) and v(t), obtained from a deterministic model - a molecular 
dynamics simulation. Note that x s and v s in ( 14.36b denote the positions and velocities 
of the particles in the probabilistic model. Then it is easy to see that if W(S; t) is 
given by 

W(S; t) = w( Xl (t) - as?) • • • w(x N (t) - x s N ), (4.37) 

then the definitions given by ( 13.131 1. (13.15b and ( 14.18b are consistent and W given by 
the above formula satisfies Liouville's equation (see ( 13.10b ). which was used in deriv- 
ing the balance equations in Section|3] Note that unlike Section[3] W(S; t) defined 
in ( 14.37b is not a probability density function. (Its integral over phase space diverges, 
since it is independent of v.) The key difference between the two approaches is that 
all quantities in the Irving-Kirkwood-Noll procedure are probabilistic, while this 
is not true for the Murdoch-Hardy procedure, if the above probabilistic interpreta- 
tion is adopted. For example, f a p in the Murdoch-Hardy procedure is deterministic. 
Therefore the structure inherent in (13.61b through the marginal densities is absent in 
the Murdoch-Hardy procedure, thus giving additional non-uniqueness. It is shown in 
Section [5] that the Hardy stress can be derived using both approaches, while the DA 
stress is a result of the Murdoch-Hardy procedure alone. 



This is only a mathematical argument. No physical significance should be drawn from this analogy. 
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4.3 Definition of the spatially-averaged traction vector 

We close this section by defining the spatially-averaged traction vector, t w (x, n; t), 
for a weighting function w, at a point x relative to a plane with normal n(x). One 
possibility is to adopt the Cauchy relation using the spatially-averaged stress tensor, 

t w (x,n;t) :— cr w (x,t)n. (4.38) 

However, since cr w is defined as a volume average, we immediately see that with 
this definition t w depends not only on the bonds that cross the surface, but also on 
nearby bonds that do not cross it J44). Hence, equation (I4.381 l does not appear to be 
consistent with Cauchy's definition of traction. 

We therefore seek an alternative definition for the spatially-averaged traction vec- 
tor. In Section l3~71 we showed that the pointwise traction vector at a point on a surface 
is the expectation of the force per unit area of all the bonds that cross the surface, mak- 
ing it a property of the surface. We would like the spatially-averaged traction to have 
the same property. We therefore define it as an average over a surface rather than 
over a volume as for the stress. For simplicity, we consider the weighting function 
Wh, defined to be constant on the averaging domain, which is taken to be a gener- 
alized cylinder of height h, with its axis parallel to n and enclosing x. The traction 
t w (x, n; t) is defined as 

t w (x,n;t) := lim tr w Jx, t)n, (4.39) 

h— ¥0 

where a Wh is the stress associated with the weighting function, Wh - In a more general 
case, an arbitrary averaging domain can be collapsed onto a surface passing through 
x, in many ways. Although this can be made mathematically more precise, we do not 
pursue that in this work. Definition ( |4.39t has a two-fold advantage over the definition 
in d438l : 

1 . The traction vector is defined to be non-local on a surface, thus making it a prop- 
erty of the surface. This is physically more meaningful, and closer to the contin- 
uum definition. 

2. The above definition differs from the traction definition in ( 14.38b . because only 
the bonds which cross the surface contribute to the traction field. 

In Section 15.21 we use definition ( 14.39b to define the Tsai traction starting with the 
spatial averaging discussed in Section 14.11 and in this way establish a link between 
the Tsai traction and the Irving-Kirkwood-Noll procedure. 



5 Derivation of different stress definitions and the issue of uniqueness 

In this section, we systematically derive various stress tensors commonly found in the 
literature from the methods developed in Section [3] and Section [4] The stress tensors 
discussed in this section are the Hardy, virial and DA stress tensors and the Tsai 
traction. 
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5.1 Hardy stress tensor 

The Murdoch-Hardy procedure described in Section[4]was independently developed 
by Murdoch ll43l and Hardy 11231 . The motivation for Hardy's study was to test the 
validity of the continuum description of phenomena in shock waves. The formulas 
suggested by Irving and Kirkwood were not useful due to the lack of knowledge re- 
garding the probability density function and the infinite series expansion in the defini- 
tion of the stress. As an alternative, Hardy used what we now term as the "Murdoch- 
Hardy procedure" to propose an instantaneous definition for stress, for the special 
case of pair potential, given by 



^( X ,t) = i £ ^W-^^^-^^^ V^^x.,^), (5.1a) 

a a " Xa X P" 

og{x,t) = -^w{x a {t) -x)m a v™\t)®v™\t), (5.1b) 

a 



where b is the bond function defined in ( 14.161 ) and is the velocity of particle a 
with respect to the continuum velocity, as defined in ( 14.241 ). To simplify the notation, 
the explicit dependence of x a and t>Jj cl on time is dropped from here onwards. Equa- 
tions ( 15. 1 ab and d5.1bl ) may look familiar. They are similar to the spatially-averaged 
generalized stress in ( 14.9b and ( 14.17b (for the special case of a pair potential). If in 
these relations, the ensemble average is replaced by a time average, we obtain a time- 
averaged Hardy stress. However, in performing such an operation, we must note the 
following: 

1. Under conditions of thermodynamic equilibrium (see footnote [8] on page ITOl). 
ensemble averages can be replaced by time averages provided that the system is 
assumed to be ergodic. Strictly speaking this time average should be done for 
infinite time, but for practical reasons we are restricted to finite time. 

2. The Hardy stress tensor is valid under non-equilibrium conditions assuming that 
the system is in local thermodynamic equilibrium^^ at all points at every instant 
of time. This is plausible only when there is a clear separation of time scales 
between the microscopic equilibration time scale r and macroscopic times. Here, 
t is not being defined rigorously. Roughly speaking, r must be sufficiently small 
so that macroscopic observables do not vary appreciably over it. 



36 Local thermodynamic equilibrium is a weaker condition than uniform thermodynamic equilibrium 
(see footnote |8). The assumption is that the microscopic domain associated with each continuum particle 
is locally in a state of uniform thermodynamic (or at least metastable) equilibrium. This is the reason why 
concepts like temperature can be defined as field variables in continuum mechanics. See for example 1171 . 
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Under these assumptions, we may replace ensemble averages with time averages in 
(l4~9l and (14.17b to obtain 

t + T 

w(x a - x)m a v™ x ® v™ l dt, (5.2a) 

t + T 

<r w , v (x,t) = — / [-f a f j ®{x a -x j3 )b{x- ) x a ,x [j )]dt, (5.2b) 

a =£f3 

where f a p, corresponding to a given potential extension, is defined in Q.511 l. and r 
represents a microscopic time scale. We see that the Hardy stress is obtained through a 
rigorous process beginning with the statistical mechanics concepts introduced in Sec- 
tion[3] From here on, we will denote the stress in (15 ,2b as the "Hardy stress", although 
we note that this definition constitutes a generalization of the original Hardy stress 
to arbitrary potentials and includes time averaging. Incidentally, the Hardy stress can 
also be derived from the systematized Murdoch-Hardy procedure described in Sec- 
tion 14.21 The kinetic part of the Hardy stress is the same as that obtained in the 
Murdoch-Hardy procedure. The potential part of the Hardy stress is derived using 
the generator g H given in ( 14.291 . This was done in Section|4](see ( 14.321 )). 

Note that the Hardy stress tensor is symmetric. One could modify this to a general 
form by choosing an arbitrary path of interaction, thus leading to a non-symmetric 
form (see Section 14.21 . Also note that the stress tensor resulting from the genera- 
tor c/ HD (see Fig. 14.3b would be symmetric, because the interaction force between 
two continuum particles is always aligned with the line connecting them. It is very 
important to observe that under non-equilibrium conditions, where we assume a lo- 
cal thermodynamic equilibrium at every instant of macroscopic time, we may assume 
that the averaging domain centered at a position x moves with the continuum velocity 
v(x,i). This fact will be used in Section IBT21 



(T Wtk (x,t) = -- y^y 



5.2 Tsai traction 

Cauchy's original definition of stress emerges from the concept of traction acting 
across the internal surfaces of a solid via the bonds that cross the surface. It is there- 
fore natural to attempt to define traction at the atomic level in a similar vein in terms 
of the force in bonds intersecting a given plane. This approach actually goes back 
to Cauchy himself as part of his effort in the 1820s to define the stress in crystalline 
systems |5]|6), which is described in detail in Note B in Love's classical book on 
the theory of elasticity [35). Cauchy's derivation is limited to zero temperature equi- 
librium where the atoms are stationary. This approach was extended by Tsai 11631 to 
the dynamical setting by also accounting for the momentum flux of atoms moving 
across the plane. The expression for the traction given in ll63l appears to be based on 
intuition. 

In this section, we show how the Tsai traction can be systematically derived from 
the Hardy stress tensor, which itself was derived from the generalized stress tensor 
defined in Section |4~TI We will see that the potential part of Tsai's original definition 
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agrees with the results of our unified framework. However, Tsai's expression for the 
kinetic part of the traction depends on the absolute velocity of the particles and there- 
fore is not invariant with respect to Galilean transformations. We show below that the 
correct expression for the Tsai traction vector t(x, n; t) across a plane P with normal 
n is 

. ( ,s _ 1 f t+T \^ , (gg -xp)-n 
t w{x ,n;t)-—J t 2^ fop \ {Xa _ Xp) . n \ dt 

J_ v m a v2\t„)(vZ\t„) ■ n) 

Ar^ p |<(M-n| ' K '> 

where r indicates the microscopic time scale, J2 a /3nP indicates the summation over 
all bonds a — f3 crossing the plane P, J^at^p indicates summation over all particles 
that cross P in the time interval [t, t + r]r I v™ 1 denotes the local relative velocity 
of particle a, and t++ indicates the time at which the particle crosses the plane. The 
correct form for v™ 1 is not immediately obvious. Below, we derive equation ( 15.3b and 
obtain an explicit expression for v™ . 

We start with the Hardy stress in (15.2b . Recall from Section l4~3l that if the averag- 
ing domain is taken to be a generalized cylinder Ch of height h, the spatially-averaged 
traction field, t w (x, n; t), on a surface passing through x with normal n is 

t w (x,n;t) = lim cr w n = lim (a w , v n + a w ,vn) (5.4) 

h— s-0 h— >0 



Using ( 15.2l i. we rewrite the potential part and kinetic part of (15.4b as 

1 f t+T 

t w y (x,n;t) = — lim / V" [-f a p ® (x a - xp)b h (x; x a , xp)]dt, (5.5a) 

IT h->0 L * — ' 
1 f t+T 

t W k(x,n;t) = lim / m a w(x a — x; h)v^(t; h)(v T ° l (t] h) ■ n)dt, 

Th ^°Jt „ 

(5.5b) 

where bh denotes the bond function for a generalized cylinder of height h. Also note 
the dependence of v™ 1 on h in ( I5.5bb . Let us first consider the potential part of the 
traction in ( I5.5ab . As h approaches zero, the generalized cylinder will no longer con- 
tain complete bonds. Assuming a constant weighting function, the bond function bh 
equals the fraction of the length of the bond lying within the generalized cylinder per 
unit volume: 

b h (x]X a ,Xp) = -^-j-. ^— r = —r- r (5.6) 

nA | (x a ~ Xfj) ■ n\ A\ (x a ~ xp) ■ n\ 



37 A particle is counted multiple times if it crosses the plane multiple times. 
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for any bond a — j3 crossing the cylinder. Therefore (I5.5ab takes the form 

ft+T 

(x a - xp) ■ n\ 



1 f t+T 
t W)V (x,n;t) = — j 



afSDP 



dt. (5.7) 



Note that the 1/2 factor is dropped because of the definition of the summation in 
the above equation. This is the first term in d5.3l ). Turning to the kinetic part of the 
traction in d5.5bb . we interchange the summation and integral to obtain 

t w k(x,n;t) = — lim } / m a w(x a — x;h)v^(t;h)(v^(t;h) ■ n)dt, 

(5.8) 

where t\ (a; h) and ti (a; h) are the times of entry and exit of particle a, respectively, 
from a cylinder of height h. The summation in the above equation is over all particles 
that are in the generalized cylinder during the time interval [t, t + r], with a particle 
counted k times if it enters and exits the cylinder k times. Multiplying and dividing 
the above equation by t-2 (a;h) — t\ (a; h) and substituting in w, we have 

J_ \- h(a; h) - h{a; h) J^ff m^jt; h)(v r a el (t; h) ■ n)dt 
wM{n > ~ At /™o 2^ h h (a; h) - Ua; h) 

= ~Tr E Ito * 8 ^^:* 1 ^ 5 ^ ^^)^^) • n), (5.9) 

where we have used the Lebesgue differentiation theorem |fl9l in the last equality. 
Note that the interchange of limit and summation in the above step is valid since we 
can assume that the summation for any Ch is a finite summation which is physically 
meaningful. Since the averaging domain moves with a continuum velocity we note 
that 

u \ , . (5.io) 

In words, this equality states that the net time spent by particle a in the cylinder, 
divided by its height, is equal to the inverse of the velocity of particle a along the 
axis of the cylinder. This is correct in the limit, h — > 0, where particles only enter and 
exit the cylinder at its ends. Substituting ( 15.10b into i5.% , we have 



At l^ cl (^) • n\ 

= ~ ™ Q < ol (^)sign« ol (^) ■ n). (5.11) 

a^P 

This is the second term in (15.3b . Note that 

v™\t) = lim v™\t; h) = v a (t) - lim v(x; h). (5.12) 

h— >0 h— >0 

Hence, we have implicitly assumed that lim/ l _ ! .o v(x; h) is well-defined for our av- 
eraging domain (plane P) which is a limit of the generalized cylinder Ch- In the 
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following calculation, we show that v(x;h) is well-defined and its exact form is de- 
rived. 

We know that for a generalized cylinder 



v(x; h) 



7 ft T J2 a m a w ( x <x(t) ~ x;h)v a (t)dt 



rt+r 
~t 2 (a;h) 



7 ft T T l p m w ( x p{ t )- x ) dt 



EJt'5) rn a v a (t)dt 



ft^-h) m P dt 
J2 a m a [x a (t 2 (a; ft)) - x a (ti(a\ h))] 
£)' ro/, [h{P;h)-h{p;h)] 



(5.13) 



where indicates summation over those particles that cross P in the time interval 
[t, t + t], including multiple entries and exits. 

Considering the limit h — >■ of the first partial fraction of the last equation we 
have 

x a (t 2 {a\h))-x a {t 1 (a-Ji)) 
jj m ° t 2 (a;/i)-t 1 (q;/i) = ^W^j „ ^ 

™° E>/» S%K|SS ~E',m,\v a (t^.n/v^).n\' 

using the fact that in the limit h -> 0, (^(Z?; ft) — ii(/3; h))/(t 2 (a; h) — h(a; ft)), 
which is the ratio of the times spent by particles f3 and a in one of their sojourns into 
the cylinder, is equal to the inverse ratio of their normal velocities. Using ( 15.14b and 
taking the limit h — > of ( 15.13b . we obtain 

v(x) = lim v(x; h) = V — m * v <*(t«) _ (5 15) 

Note that the above expression for the continuum velocity is far from intuitive. One 
might expect the continuum velocity to be the average velocity of particles crossing 
the surface, but this is not true. It is clear from the above equation that the averaging 
is not trivial. 

From the relationship between the Tsai traction in (15.3b and the Hardy stress 
tensor in d5.2l ), it is apparent that the Tsai traction is a more local quantity than the 
Hardy stress tensor. The Tsai traction performs better than the Hardy stress in systems 
with free surfaces. This was studied by Cheung and Yip Q for a one-dimensional 
case, in which virial stress and Tsai stress are compared (the virial stress is a special 
case of Hardy stress as shown in the next section). 

The Tsai traction definition can be used to evaluate the stress tensor at a point by 
evaluating the traction on three perpendicular planes(3 However, it is not clear from 
the perspective put forward by Tsai 11631 whether the resulting stress tensor would 



38 For example, if the normals to the planes are aligned with the axes of a Cartesian coordinate sys- 
tem with basis vectors e;, then t(ei) would give the components an, 021, &31, t(e2) would give the 
components <ri2, o"22, C32, and t(es) would give the components 013, <J23, C33. 
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be symmetric or even well-defined, i.e., it is not clear if another choice of planes 
will give suitably transformed components of the same stress tensor. Our derivation 
suggests that a stress tensor constructed from the Tsai traction should be well-defined 
and symmetric, at least in a weak sense, since it is a limit of the Hardy stress, which 
has these properties. The numerical experiments presented in Section[6l suggest that 
the Tsai traction is invariant with respect to the position of the Tsai plane P and the 
resulting stress tensor is symmetric. 



5.3 Virial stress tensor 

In this section, we show that the virial stress tensor derived in Section [2] and in Ap- 
pendix[A]can be re-derived from the time-averaged version of the Hardy stress given 
in ( 15.2b . The expression for the virial stress tensor is obtained from ( 15.2b as a special 
case for a weighting function which is constant on its support. The bond function, b, 
in (15.2b is evaluated approximately using its definition (14.16b by only counting those 
bonds a — (3 that lie entirely within the averaging domain and neglecting the bonds 
that cross the averaging domain. Hence, b(x; x a , xp) is given by 

where f2 x denotes the averaging domain centered at x. Substituting ( 15.161 ) into ( 15.21 ), 
we have 



1 f t+T 1 

cr ( 3; ^) = TVQ u n ) / ~ E m «< el ®< el +2 E [-/a0®OBa-SB0)] 

(5.17) 

which is identical to ( IA.27b in Appendix lAl It is clear from this that the virial stress 
tensor is only an approximation and tends to the Hardy stress as the volume of the 
averaging domain is increased. This is because the ratio of the measure of bonds that 
cross the surface to those which are inside the averaging domain decreases as the 
size of the domain increases. The difference between the virial stress tensor and the 
Tsai traction was analytically calculated for a one-dimensional chain by Tsai (see 
ll63l ). Since, taking the averaging domain size to infinity is equivalent to taking the 
thermodynamic limit in this context, the Hardy and virial stress expressions become 
identical in this limit. Since the virial theorem was also derived in Section [2] for the 
case of equilibrium statistical mechanics, it follows that the Irving-Kirkwood-Noll 
procedure is consistent with the results of equilibrium statistical mechanics in the 
thermodynamic limit. 



dt. 



5.4 DA stress tensor 

It was seen in Section l4~2l that the DA stress tensor, defined in ( 14.33b . is derived using 
an appropriate generator in the systematized Murdoch-Hardy procedure. However, 
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unlike the Hardy stress, the DA stress cannot be derived from the Irving-Kirkwood- 
Noll procedure. It is also worth noting that the stress tensor given by ( 14.33b is in 
general non-symmetric, and only under very special conditions yields a symmetric 
tensor B31l . 



5.5 Uniqueness of the macroscopic stress tensor 

Three possible sources of non-uniqueness for the stress tensor have been identified in 
our discussion: 

1. Given that there are multiple potential extensions (see Page l24li. different force 
decompositions are possible and hence different pointwise stress tensors can be 
obtained. 

2. For a given pointwise stress tensor, a new pointwise stress, which also satisfies 
the balance of linear momentum, can be obtained by adding on an arbitrary tensor 
field with zero divergence. 

3. The generalization of the Irving-Kirkwood-Noll procedure in Section [3~6l to arbi- 
trary "paths of interaction" leads to the possibility of non-symmetric expressions 
for the pointwise stress tensor. 

We address the first two issues in this section. The third source of non-uniqueness 
is only possible in systems where the discrete particles making up the system pos- 
sess internal structure, such as internal polarization or spin. For systems of discrete 
particles without internal structure only straight bonds are possible due to symmetry 
arguments. We leave the discussion of particles with internal structure to future work. 



Uniqueness and potential energy extensions 

The first source of non-uniqueness of the stress tensor is related to the potential en- 
ergy extension discussed in Section [3~4l We show below that the macroscopic stress 
tensor, calculated as a spatial average of the pointwise stress tensor with constant 
weighting function, is always unique in the thermodynamic limit (see footnote Q] on 
page |5]l, i.e., the difference between the spatially-averaged pointwise stress tensors 
resulting from two different extensions tends to zero, as the volume of the averaging 
domain is increased. 

The discussion below is limited to 5-body potentials since it can be easily ex- 
tended to any interatomic potential. We first show that the contribution due to any 
cluster of 5 particles within the averaging domain is zero. Without loss of generality, 
we may assume that our system consists of 5 particles interacting with an interatomic 
potential energy given by 



V int = V(a;i,...,a;5). 



(5.18) 
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Fig. 5.1 A cluster of 5 particles that lie completely inside the averaging domain, does not contribute to the 
ambiguity in the stress tensor. 



Let Vint(Ci2i ■ • • j C45) an d Knt(Ci2, • • • , Cis) t> e tw0 different extensions of Vint from 
the shape space S to R 10 (see Section [3~4l i, and for any s = (ri2, . . . , ^45) 6 S, let 

f a p(x 1 ,...,x s ) := —— (s)— , (5.19) 

f* fi{ x 1 ,...,x 5 )^^(s)^^, (5.20) 

P OQa/3 r a p 

be their corresponding force decompositions. Let cr and cr* denote the resulting 
pointwise stress tensors in the Irving-Kirkwood-Noll procedure from Vi nt and V* nt , 
respectively. Let fl x denote the averaging domain^ centered at x that is used to cal- 
culate the Hardy stress tensor. Using (I5.2bb and noting that all the bonds lie within i?, 
the difference between the Hardy stress tensors resulting from these two representa- 
tions, for the special case of a constant weighting function, is given by 

1 f t+T 

Acr(x, t) := (T W -(T* W = — — ^2 / \- A fa& ® ( x a - xp)]dt, (5.21) 

where Af a0 := f a/3 - f* p . 

We would like to show that A<jri\ = 0, where ri\ is the normal vector as shown 
in Fig. 15.11 The essential idea to is to interchange the integration and summation in 
J5.2U and split the terms appearing in the summation into fractions, such that each 
fraction yields a zero contribution to Acrn\ . In order to show this, we partition the av- 
eraging domain into regions such that no region contains a particle in its interior and 
the partition surfaces are perpendicular to the normal (see Fig. 15.11 ). Let hp denote 
the width of the partition P. Using (15.211) . we can now write Acrni as 



t+T _ 

-Af al3 (x a - xp) ■ m 1 



p a/3np 

t+T 



2r vo 



x a - x f3 ) ■ ni| 



For simplicity assume that the averaging domain is convex. 
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where J2 a pnP denotes the summation over the bonds crossing partition P, and 

(x Q - xp) ■ ni 



AF P 



E 

Q /3nP L 



x a - xp) ■ txi| 



(5.23) 



is the net force on particles on one side of the partition due to particles on the other 
side. Since both representations give the same total force on each particle, the force 
difference, or net force, on each particle is zero and therefore, AFp = 0. For exam- 
ple, for the partition shown in the figure, 



AF P 

Since 4/45 = — zA/54, we have 



-2(4/ 5 i+4/ 43 ). 



(5.24) 



AF P = -2(Afsi + 4/43 + 
= -2(^/43 + 4/45) - 
= -24/f* - 2Af 5 nt 



4/45 + 4/ 54 ) 
- 2(Z\/ 5 i + 4/54) 
= + = 0. 



(5.25) 



Hence, Acrni = 0. Undertaking a similar argument in the other directions, we see 
that Acrni = 0. These results together imply that Act = 0. Given this, we can con- 
clude that any cluster of particles that lies entirely within the averaging domain does 
not contribute to the spatial average of the difference between two stress definitions. 
Consequently, the only non-zero contribution comes from those clusters for which 
the bonds connecting its particles cross the averaging domain. Since this contribution 
scales as surface area, it tends to zero as volume tends to infinity. 



Uniqueness and the addition of a divergence-free field to the stress 

The second source of non-uniqueness of the stress tensor involves the addition to it 
of a divergence-free field. This issue is partly addressed by the result (shown in Sec- 
tion [53} that the spatially-averaged pointwise stress converges to the virial stress in 
the thermodynamic limit (see footnoteQ]on page[5}. Consider the pointwise stress, cr, 
obtained through the Irving-Kirkwood-Noll procedure, which satisfies the balance 
of linear momentum, and a new pointwise stress, cr = cr + cr, where div x cr = 0. 
Clearly, cr also satisfies the balance of linear momentum and is therefore also a valid 
solution. The spatially-averaged stress obtained from the new definition is 

& w (x,t)= w(y-x)cr(y,t)dy = w(y-x)(a(y, t)+a(y, t)) dy. (5.26) 

We showed in Section 15.31 that in the thermodynamic limit, the spatially-averaged 
pointwise stress, cr, converges to the virial stress. We also expect & w to equal the 
virial stress in this limit (since any macroscopic stress must converge to this value 
under equilibrium conditions). Therefore, (15.26b reduces to 

lim / w(y — x)&(y, t) dy = 0, (5.27) 

where limxD refers to the thermodynamic limit. Equation ( 15.27b places a strong con- 
straint on allowable forms for cr, the implications of which are left for future work. 
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Fig. 6.1 The virial pressure as a function of time is plotted for an isolated cube of aluminum at 300K. The 
total pressure/virial pressure is the sum of the kinetic and potential pressures. 



6 Numerical Experiments 

In this section, we describe several numerical experiments, involving molecular dy- 
namics and lattice statics simulations, conducted to capture differences in the spatially- 
averaged stress measures derived in Section [5] We consider the Hardy stress defined 
in ( 15.2b . the Tsai traction defined in ( 15.3b . the virial stress defined in ( 15.171 ) and the 
DA stress defined in ( 14.331 1. We will sometimes refer to these as the "microscopic 
definitions" or the "microscopically-based stress tensors". 



6.1 Experiment 1 

We begin with the study of the kinetic part of the stress tensor. From the discussion 
in Section [5] it is clear that unlike the definition for the potential part of the stress 
tensor, there is no ambiguity in the definition for the kinetic part of stress. However, 
the kinetic part of the stress may appear to be at odds with the continuum definition 
of stress that is stated solely in terms of the forces acting between different parts of 
the body. The need for the kinetic part of stress becomes apparent when considering 
an ideal gas, where the potential interaction term is zero by definition and therefore 
it is the kinetic term that is wholly responsible for the transmission of pressure. 

To demonstrate that the kinetic term in the stress tensor does indeed exist, we 
perform the following constant energy molecular dynamics simulation of an isolated 
cube. The cube, consisting of 4000 aluminum atoms in a face-centered cubic (fee) 
arrangement (10 x 10 x 10 unit cells), is floating freely in a vacuum. The atoms interact 
according to an EAM potential for aluminum due to Ercolessi and Adams ifTBI : 

a,P a p 

a^P P^a 
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Here U a , called the embedding function, is the energy required to embed particle a in 
the electron density, p a , due to the surrounding particles, and fp{r a p) is the electron 
density of particle j3 at x a . The initial positions of the atoms are randomly perturbed 
by a small amount relative to their zero temperature equilibrium positions and the 
system is evolved by integrating the equations of motion. The initial perturbation 
is adjusted so that the temperature of the cube is about 300K (small fluctuations in 
temperature are expected since temperature is not controlled in the simulation). Since 
the block is unconstrained, we expect the stress, er, in the box and consequently the 
pressure, defined by p = tr er, to be zero. The virial expression for calculating 
the pressure follows from (IA.231 > as 



P = 



1 

W 



E 



|E1 

a^0 



fa 



(6.2) 



where 



fa/3 



EAM 



Or, 



a0 



(6.3) 



The three curves shown in Fig. 16. H are the potential and kinetic parts of the pres- 
sure and the total pressure as a function of time, calculated using ( 16.2b . As expected 
the total pressure tends to zero as the system equilibrates. However, the potential 
and kinetic parts are non-zero, converging to values that are equal and opposite such 
that their sum is zero. More interestingly, the kinetic part is not insignificant for our 
system. This clearly shows that kinetic part cannot be neglected even when consid- 
ering solid systems. This can be quantified by noting that the kinetic part in ( 16.21 ) 
is simply the temperature per unit volume given by the equipartition theorem |26i], 
ksT = 2T/3N, where T is the kinetic energy. Therefore (16.21 i reduces to 



1 

P= V 



Nk 



1 R a 



(6.4) 



For example at 300 K, k^T = 0.02585 eV. The lattice spacing for the system consid- 
ered is equal to 4.032A. Hence, the volume per atom is V/N = 4.023 3 /4 = 16.387A 
and the kinetic pressure is 1.577 meV/A. This translates to 252.394 MPa, which is a 
considerable stress. 



6.2 Experiment 2 

It is clear from Section [64] that the kinetic stress is a sizable quantity and cannot be 
neglected. In this experiment, we further explore the interplay between the potential 
and kinetic parts of the stress. 

Consider a crystalline solid at a relatively low temperature under uniform stress. 
The atoms will vibrate about their mean positions with an amplitude that is small rel- 
ative to the nearest-neighbor spacing. Now imagine placing a Tsai plane P between 
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Fig. 6.2 The effect of the position of the Tsai plane on the potential and kinetic parts of stress. Frames (a) 
and (b) show schematic diagrams of a two-dimensional triangular lattice with (a) the Tsai plane positioned 
midway between the lattice planes and (b) the Tsai plane positioned almost on top of a lattice plane. The 
open circles correspond to the ideal lattice positions. The black circles are atoms that are shown in mid- 
vibration relative to the lattice site as indicated by the arrow. The Tsai plane is indicated by a vertical 
dashed line. The bonds crossing it appear as dotted lines. Frame (c) shows the plot of the kinetic part of 
stress Oj,, potential part of stress crj 1 and the total stress an, as a function of the normalized position 
sp = (xp —xl) I Ax of the Tsai plane P, where xp is the position of P, xl is the position of the lattice 
plane, and Ax is the spacing between the lattice planes. 



two crystal lattice planes and measuring the traction across it. If P is midway between 
the lattice planes (see Fig. |6.2(a)| i, we expect that relatively few atoms will cross P 
and that consequently the kinetic stress will be small or even zero. In contrast, if P 
is close to the lattice plane there will be many such crossings and the kinetic stress 
will be large in magnitude. This seems to suggest that the traction will change as a 
function of the position of P, which would be incorrect since the system is under uni- 
form stress. The reason that this does not occur is that every time an atom crosses P, 
the bonds connected with it reverse directions with respect to P, changing a positive 
contribution to the contact stress to a negative one and vice versa (see the bonds con- 
nected with atom A in Fig. |6.2(a)| and Fig. |6.2(b)| i. This effect on the potential part 
of the stress exactly compensates for the change in magnitude of the kinetic stress 
leaving the total stress constant. This is demonstrated numerically in Fig. |6.2(c")| This 
graph shows the results obtained from a molecular dynamics simulation of the sys- 
tem described in Section l6"Tl with periodic boundary conditions. The periodic length 
of the box is set based on the zero temperature lattice spacing. Consequently upon 
heating by a temperature change of AT, a compressive stress is built up in the box 
according to 

e = s : er + 1 cut AT = 0, (6.5) 

where s is the elastic compliance tensor and ar is the coefficient of thermal expan- 
sion. Inverting this relation for an fee crystal with cubic symmetry oriented along the 
crystallographic axes, we have 



oil = 0-22 = 0-33 = -(en + 2c 12 )AT = a, 



(6.6) 
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with the rest of the stress components zero. In ( 16.6b . Cy are the elastic constants of the 
material. Substituting in the appropriate values for Ercolessi-Adams EAM aluminum 
fl6l (en = 118.1 GPa, ci2 = 62.3 GPa, a T = 1.6 x lCT 5 !^ 1 ) and Z\T = 310K, 
gives a = —1.2 GPa. We see that the total stress in Fig. |6.2(c)| is constant regardless of 
the position of the Tsai plane and equal to the expected value of — 1 .2 GPa computed 
above. However, the kinetic and potential parts change dramatically. When the Tsai 
plane is away from the lattice planes (sp = ±0.1), the kinetic stress is zero and the 
entire stress is given by the potential part of the stress. As the Tsai plane moves closer 
to a lattice plane (\sp\ —> 0), the kinetic stress becomes more negative (increasing 
in magnitude) and the potential part of stress increases in exactly the right amount 
to compensate for this effect. When the Tsai plane is right on top of a lattice plane 
(sp = 0), both the kinetic stress and potential stress are maximum in magnitude, but 
their sum remains equal to the constant total stress. This is a striking demonstration 
of the interplay between the kinetic and potential parts of the stress tensor. 



6.3 Experiment 3 and 4 

In this section, the predictions of the microscopically-based stress tensors are com- 
pared with analytical solutions from elasticity theory for two simple boundary-value 
problems. This is a revealing test, since stress is a continuum concept and therefore 
the microscopic definitions should reproduce the results of a continuum theory under 
the same conditions. We perform two numerical experiments. In each experiment, an 
atomistic boundary-value problem is set up, and the values computed from the dis- 
crete system are compared with the "exact" result computed from elasticity theory 
for the same problem using material properties predicted by the interatomic potential 
used in the atomistic calculations. The numerical experiments are conducted at zero- 
temperature since there is no controversy regarding the form of of the kinetic stress 
which is the same for all stress definitions. Therefore, a comparison at zero tempera- 
ture is sufficient to probe the differences between the stress measures, at least under 
equilibrium conditions. The properties we are interested in studying are: 

1 . Symmetry of the stress tensor. 

2. Convergence of the stress tensor to the continuum value with the size of the av- 
eraging domain (a three-dimensional volume in the case of virial, Hardy and DA 
stresses and a plane in the case of the Tsai traction). 



Inter-atomic Model 



The numerical experiments in this section are carried out using a Lennard-Jones po- 
tential. The exact choice of material parameters is unimportant, since the objective 
of the experiment is to compare the values obtained from the microscopically-based 
stress for the discrete system with the "exact" values obtained from the continuum 
elasticity theory for the same material. The Lennard-Jones parameters, e and a, are 
therefore arbitrarily set to 1 . The potential has the following form: 



1 



- 0.0078^ + 0.0651. 



(6.7) 
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Note that the above equation has been rendered dimensionless by expressing lengths 
in units of a and energy in units of e. As seen in the above equation, the Lennard- 
Jones potential is modified by the addition of a quadratic term. This is done to ensure 
that </>(r cu t) = and </>'(r cut ) = 0, where r cut = 2.5, denotes the cutoff radius 
for the potential. We refer to this as the "modified Lennard-Jones potential". The 
ground state of this potential is an fee crystal with a lattice constant of a = 1.556 
and elastic constants, en = 87.652, c\i = C44 = 50.379. The conventional elastic 
moduli associated with the cubic elastic constants are 041 : 

E = (c 2 n + c u c 12 - 2c? 2 )/(c u + ci 2 ) = 50.877, (6.8) 
,u = c 44 = 50.379, (6.9) 
v = c 12 /(c n + C12) = 0.365, (6.10) 

where E is Young's modulus, /1 is the shear modulus, and v is Poisson's ratio. (In the 
above, elastic constants are given in units of e/er 3 . Poisson's ratio is dimensionless.) 

Experiment 3: Dependence of the microscopically-based stress on the averaging do- 
main size 

The main aim of this experiment is to study the dependence of the stress given by var- 
ious definitions (Hardy, Tsai, virial and DA) on the size of the averaging domain. We 
consider the special case of uniform uniaxial loading with an = 1 (all other stress 
components zero). Our system is a cube of 10 x 10 x 10 unit cells (4000 atoms) with 
periodic boundary conditions applied in all directions. To impose the uniaxial load- 
ing, the periodic lengths li (i = 1, 2, 3) in the three directions are modified according 
to the linear elastic solution for uniform straining: 

h = l0a(l + a n /E) = 10.197, (6.11) 
h = h = 10a(l - vaix/E) = 9.928. (6.12) 

We then compute the stress at the center of the periodic cell, while increasing the size 
of the averaging domain. In comparing the different stress definitions, the domain 
size is set by the Tsai plane which is taken to be a square normal to the 1 -direction 
with the same dimension w in the 2 and 3-directions. The averaging domain for the 
Hardy, virial and DA stresses is a sphere of diameter d = w. The weighting function, 
w(r), for the Hardy stress is taken to be constant with a suitable mollifying function, 

!c if r < R — e 

\c [1 - cos (^tt)] if R- e < r < R , (6.13) 
otherwise 

where c is chosen appropriately to normalize w. The results are presented in Fig. 16.31 
where the stress 011 is plotted as a function of the normalized domain size, s — 
w/10a. (Recall that the applied value is <th = 1.) We make the following qualitative 
observations based on these results: 

1 . The Hardy stress converges to the exact value most quickly of all stress definitions 
and has the least noise. 
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Fig. 6.3 Plot showing the dependence of crn, calculated using different definitions on the averaging do- 
main size. The variable s represents the ratio of the domain size to the length of the system (10 unit cells). 

2. The normal stress computed from the Tsai traction oscillates about the exact value 
with a fluctuation amplitude that decays rather slowly with domain size. The os- 
cillations reflect the symmetry of the crystal as new bonds enter the calculation 
with increasing plane size. 

3. The virial stress is always smaller than the Hardy and Tsai stresses since it does 
not take into account the bonds that cross out of the averaging domain. It appears 
to be converging towards the exact value, but convergence is slow and even at the 
maximum domain size studied, the virial stress still has a significant error. 

4. The DA stress is much smaller than all other stresses due to greater averaging. 

Experiment 4: A plate with a hole under tension 

We now consider one of the classical elasticity boundary-value problems: an infinite 
plate with a hole subjected to uniaxial tension at infinity. This is traditionally 
named the Kirsch problem for an isotropic material model. Our objective is to com- 
pare the microscopically-based stresses computed for a discrete system set up for the 
Kirsch problem with the exact solution. A complication in making this comparison 
is that the fee Lennard- Jones material we are considering is crystalline with cubic 
symmetry and is not isotropic. We must therefore compare the discrete solution with 
the more general solution for the Kirsch problem from the theory of elasticity for 
anisotropic media 041 . For anisotropic materials, the stress concentratior@ at the 
hole is no longer 3 (as it is for an isotropic material), but depends on the elastic con- 
stants of the material. For the elastic constants of the Lennard- Jones model in (16.7b . 
we obtain a stress concentration of 2.408. In addition to the overall stress concentra- 
tion, the analytical solution provides the complete stress field about the hole. We can 
therefore compare the microscopically-based stress fields with the continuum result. 



40 The stress concentration is defined as the ratio of the maximum stress to the applied stress (Too ■ The 
maximum stress for the Kirsch problem occurs at the circumference of the hole. 
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Fig. 6.4 Normalized an component of the stress along the x\ = line for an anisotropic plate with 
a hole subjected to uniaxial tension in the 1-direction. The X2 coordinate is normalized by the height h 
of the atomistic model. The exact solution for an infinite plate obtained from anisotropic linear elasticity 
(black solid line) is compared with the results obtained from the three microscopic definitions in the three 
columns: Hardy, Tsai and virial. The four rows correspond to four different averaging domains constituting 
1%, 2.5%, 5% and 10% of h. 



In order to model an infinite elastic space, we consider a large square plate ori- 
ented along the crystallographic axes consisting of 367,590 atoms, with a hole of 
radius 25a, where a is the lattice constant. The plate is constructed by stacking 
100 x 100 x 10 unit cells and excluding the atoms that lie withing the radius of the 
hole. The relatively large system size helps to ensure that the variation of the contin- 
uum stress is small on the lengthscale of the lattice spacing and minimizes boundary 
effects near the hole. The atoms interact according to the modified Lennard- Jones 
potential given in ( 16.71 ). 

As before, the averaging domain size is set by the length and width of the Tsai 
plane, with the Hardy, virial and DA stress using a sphere of diameter equal to the 
length of the Tsai plane. The system is loaded by an = by displacing the 
atoms according to the exact solution from continuum mechanics for linear elastic 
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anisotropic media 11341 . The applied stress is sufficiently small, so that the assump- 
tion of material linearity and the small strain approximation inherent in the elasticity 
theory provide a good approximation for the behavior of the system. After the atoms 
are displaced, the stress tensor is evaluated on a uniformly-distributed grid of points 
on the mid-plane of the plate located at X3 = 0. A grid of 100 x 100 points is chosen 
to evaluate the virial, Tsai and Hardy expressions and a grid of 30 x 30 is chosen 
to evaluate the DA stress tensor. A coarser grid is used for the DA stress due to the 
higher computational cost of this calculation. 

First, we consider the 011 component of the stress along the x\ = line, where 
we expect the maximum value at the hole surface. The results are plotted in Fig. 16.41 
which shows a comparison between the exact value and the three microscopic stress 
definitions (Hardy, Tsai and virial) for four different averaging domains ranging from 
1% to 10% of the height h of the atomistic model. We see that the Hardy and Tsai 
stresses faithfully follow the exact solution, but then drop-off as their averaging do- 
main overlaps with the hole. This drop-off reflects the fact that the microscopically- 
based stress measures are bulk expressions. The smaller the averaging domain, the 
closer the microscopic measures can approach the exact stress concentration at the 
hole surface, however, this increased fidelity comes at the cost of significantly large 
fluctuation about the exact value. The virial stress is identically zero for the smallest 
averaging domain because it is too small to contain complete bonds. For the same 
reason, the Hardy stress experiences very large fluctuations and a nearly constant av- 
erage value. For larger averaging domains, the Hardy stress has smaller fluctuations 
than the other stress definitions. 

The reason that the drop-off effect described above is so pronounced in this sim- 
ulation, is that the system is very small by continuum standards. If instead of a hole 
with a radius of 25a, we studied a plate with a hole 100 or 1000 times larger, using 
the same-sized averaging domain, the spatially-averaged expressions would get much 
closer to the correct value before dropping off over the same lengthscale as seen in 
Fig. 16.41 However, microscopic stress measures are often computed for small sys- 
tem sizes and therefore the difficulties presented in the figure are typical of realistic 
atomistic simulation. 

Next, we explore the stress field over the entire plane. The color density plots 
given in Fig. 16.51 show variation of an in the mid-plane of the plate. It can be seen 
that the stress within the hole is zero. Comparing the plots for o\\ of the virial stress 
(Fig. 6.5(b)| i, Tsai stress (Fig. |6.5(c)| i, Hardy stress (Fig. |6.5(d)| i and the DA stress 
(Fig. 6.5(e) 1 with the exact solution (Fig. |6.5(a)) , we see that the first three defini- 
tions capture the overall variation in an, whereas the DA stress does not. However, 
it is clear that the microscopically-based stress in all of the cases is smeared relative 
to the stress given by the exact solution and none reach the exact stress concentration 
of 2.408. This is a result of the averaging procedure involved in all the definitions 



as explained above. Although the DA stress tensor plotted in Fig. |6.5(e j^] captures 



the variations in the field, it is much smaller in magnitude compared to the exact so- 
lution. This is because of the greater degree of averaging involved in the DA stress 



41 Fig. |675(e)| and Fig. l6.8l are generated from a much coarser grid compared to the other plots due to 
the computational expense of the DA stress definition. 
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(b) (c) (d) (e) 

Fig. 6.5 Color density plots of en are plotted on a common scale: (a) exact (b) virial (c) Tsai (d) Hardy 
(e) DA. Results plotted for an averaging domain size of 10% of the height of the model. 




(a) (b) (c) 



Fig. 6.6 Color density plots of error in an, defined as the absolute value of (an — <r™ act )/cr°^ act , 
where an is the stress calculated using (a) virial (b) Tsai and (c) Hardy stress definitions. 

tensor. Overall, the Hardy stress is less noisy than the virial or Tsai definitions due to 
the smoothing afforded by the weighting function. (This is hard to see in the figure.) 
The stress computed from the Tsai traction, in particular, is more noisy since the av- 
eraging is limited to a plane compared with the volume averaging of the Hardy and 
virial definitions. However, this more localized definition enables the Tsai stress to 
approach the exact stress concentration most closely of all of the microscopic defini- 
tions. Similar results were observed by Cheung and Yip J7) for the stress near a free 
surface. 

The relative error in an for the three microscopic definitions is shown in Fig. 16.61 
Of the three definitions, the stress computed from the Tsai traction is generally more 
accurate, followed by the Hardy stress and then the virial stress. As noted above, the 
Tsai stress does particularly well in capturing the variations in the stress field close 
to the hole where the fact that it is localized in one direction is particularly helpful. 
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It is also interesting to examine the shear stress components. Fig. 16.71 shows the 
exact result from the continuum solution and the o\i and 021 components of the stress 
tensor computed from the Tsai traction from two different planes, one normal to the 
1 direction and the other normal to the 2 direction. We see that the Tsai stress repro- 
duces the exact distribution and appears generally symmetric. This suggests that the 
symmetry of the Hardy stress is preserved while taking the limit to arrive at the Tsai 
traction. The DA stress tensor is in general non-symmetric [45], but from Fig. |6.8(a)| 
and Fig. |6. 8(b)] we observe that in this case, an and a%y appear similar. Interestingly, 
in contrast to the normal stress, the magnitude of shear stress is captured by the DA 
stress definition, at least for the case studied here. The reason for this is not obvious. 

Overall, we can summarize our results as follows. Of the three definitions studied, 
the Hardy stress is generally preferred. It tends to be the smoothest and provides good 
accuracy away from surfaces as long as the lengthscale over which the continuum 
fields vary is large relative to the atomic spacing. In situations where either of those 
conditions break down, the Tsai traction provides a better localized measure of stress. 
The virial stress is less accurate than both. From a computational standpoint, the 
virial stress has the advantage of being easiest to compute. The evaluation of the 
bond function in the Hardy stress makes it slightly more expensive to compute, but 
comparable to the virial stress. The Tsai traction is most difficult and time consuming 
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to compute, since it requires the detection of bonds and atoms that cross a given plane 
during the averaging process. Furthermore, this evaluation must be performed for 
three separate planes in order to obtain the full stress tensor in three dimensions. 



7 Summary and future work 

In this paper, we provide a unified interpretation and possible generalization of all 
commonly used definitions for the stress tensor for a discrete system of particles. The 
macroscopic stress in a system under conditions of thermodynamic equilibrium are 
derived using the ideas of canonical transformations within the framework of classical 
statistical mechanics. The stress in non-equilibrium systems is obtained in a two-step 
procedure: 

1. The Irving-Kirkwood-Noll procedure |27 4<S| is applied to obtain a pointwise 
(microscopic) stress tensor. The stress consists of a kinetic part cry and a poten- 
tial part er v . The potential part of the stress is obtained for multi-body potentials 
which have a continuously differentiable extension from the shape space of the 
system to a higher-dimensional Euclidean space@ This generalizes the original 
Irving-Kirkwood-Noll approach that was limited to pair potentials. This general- 
ization is obtained based on the important result that for any multi-body potential 
with a continuously differentiable extension, the force on a particle in a discrete 
system can always be expressed as a sum of central forces. In other words, the 
strong law of action and reaction is always satisfied. 

2. The pointwise stress obtained in the first step is spatially averaged to obtain the 
macroscopic stress. 

This two-step procedure provides a general unified framework from which various 
stress definitions can be derived including all of the main definitions commonly used 
in practice. In particular, it is shown that the two-step procedure leads directly to the 
stress tensor derived by Hardy in ll23l . The traction of Cauchy and Tsai [5 6]|63) is 
obtained from the Hardy stress in the limit that the averaging domain is collapsed to 
a plane. The virial stress of Clausius and Maxwell [8 ,40,41 1 is an approximation of 
the Hardy stress tensor for a uniform weighting function where bonds that cross the 
averaging domain are neglected. The Hardy stress and virial stress become identical 
in the thermodynamic limit. In this manner, clear connections are established between 
all of the major stress definitions in use today. 

The unified framework described above yields a symmetric stress tensor for all 
interatomic potentials which have an extension, when used with the standard Irving- 
Kirkwood-Noll procedure. However, there are materials in nature, such as liquid 
crystals, which can have non-symmetric stress tensors. In order to explore the pos- 
sibility of non-symmetric stress, the Irving-Kirkwood-Noll procedure is generalized 
to curved paths of interaction as suggested in J52). This involves the generalization 
of Noll's lemmas in 11481 . originally derived for straight bonds, to arbitrary curved 
paths as defined in this paper. These generalized lemmas, lead to a non-symmetric 



Most practical interatomic potentials satisfy this conditions. 
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stress tensor when applied within the Irving-Kirkwood-Noll procedure. It is pos- 
tulated that curved paths of interaction may be important in systems with internal 
degrees of freedom, such as liquid crystals and objective structures ll28ll . This is left 
for future work. 

One of the key points addressed in this paper is the uniqueness of the stress ten- 
sor. Three possible sources of non-uniqueness of the stress tensor are identified and 
addressed: 

1. Different pointwise stress definitions can be obtained for different potential en- 
ergy extensions. This is demonstrated through a simple one-dimensional exam- 
ple. We also show that regardless of the uniqueness of the pointwise stress tensor, 
the macroscopic stress tensor obtained through a procedure of spatial averaging 
is unique since the difference resulting from alternative pointwise stress tensors 
tends to zero as the volume of the averaging domain is increased. 

2. The pointwise stress tensor is obtained by solving the balance of linear momen- 
tum, div x cr v = h(x), where er v is the potential part of the stress tensor, and h(x) 
is a known function. The Irving-Kirkwood-Noll procedure leads to a closed-form 
solution to this problem. However, an arbitrary tensor field <x with zero divergence 
can be added to er v without violating the balance of linear momentum. We argue 
that in the thermodynamic limit, the non-equilibrium stress obtained through our 
unified two-step process must converge to the virial stress of equilibrium statis- 
tical mechanics. This is similar to the argument made by Wajnryb et al. Il64l . 
This condition is satisfied by the general stress expression that we obtain. Any 
divergence-free stress er added to this stress must therefore also disappear under 
equilibrium conditions. This greatly restricts the allowable forms of er. 

3. The generalization of the Irving-Kirkwood-Noll procedure from straight bonds 
to arbitrary curved paths of interaction implies the existence of multiple stress 
tensors for a given system. However, the existence of curved bonds implies the 
existence of internal structure for the discrete particles, a possibility already dis- 
cussed by Kirkwood in ll29l . For a system of point masses without internal struc- 
ture, only straight bonds are possible due to symmetry arguments, and therefore 
this source of non-uniqueness is removed. The general case of non-symmetric 
stress must be addressed within the context of an appropriate multipolar theory 
as discussed by Pitteri |49l . We leave this to future work. 

In addition to the unified framework described above which is based on the 
Irving-Kirkwood-Noll procedure, we also investigated the Murdoch-Hardy proce- 
dure 112311431 of defining continuum fields as direct spatial averages of microscopic 
quantities. We demonstrate that this approach can be systematized by adopting a 
non-local continuum perspective and introducing suitable generator functions. The 
various stress definitions resulting from the Murdoch-Hardy procedure, such as the 
Hardy, virial and the "double-averaged" stress (suggested by Murdoch in j47l ) can 
be derived from this unified framework. Although we share the concern regarding 
the ambiguity of the probability density functions used in the Irving-Kirkwood-Noll 
procedure that led Murdoch to develop the direct spatial averaging approach \A6\, 
we feel that since these probability density functions exist in principle, the Irving- 
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Kirkwood-Noll formalism is the correct framework to phrase the problem in with 
approximations introduced later to derive practical expressions. 

Finally, numerical experiments involving molecular dynamics and lattice statics 
simulations are conducted to study the various stress definitions derived in this paper. 
It is generally observed that the Hardy stress definition appears to be most accurate 
and converges most quickly with the averaging domain size. In situations where a 
more localized measure of stress is needed, such as near surfaces or defects, the Tsai 
traction can be used instead. The virial stress is less accurate than the other two def- 
initions and converges most slowly with averaging domain size. Its main advantage 
is its simple form and low computational cost. One of the most interesting results, 
which requires further study, comes from Experiment 2 of a crystalline system under 
uniform hydrostatic stress. Fig. |6.2(c")| shows that although the potential and kinetic 
parts of the Tsai traction largely depend on the position of the Tsai plane between two 
adjacent lattice planes, the total stress remains constant. This calculation provides a 
striking demonstration of the interplay between the kinetic and potential parts of the 
stress tensor. 
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A Derivation of the virial stress tensor 

The original virial theorem was a scalar equation credited to Clausius Q5], which gives a definition for the 
pressure in a gas. Maxwell 1401 extended this result to a tensor version which gives a definition for stress. 
We present here a more modern version of the virial theorem partly based on the derivation by Marc and 
McMillan (38j. 

Consider a system of N interacting point masses. The position of each point mass is given by 

x a = x + r a , (A.l) 

where x is the position of the center of mass of the system of particles and r a is the position of each point 
mass relative to the center of mass. From Newton's second law, we have 

fa=Pa, (A.2) 

where f a is the force acting on particle a, and p a is its linear momentum given by 

Pa = rn a (x + r a ) = m a (x + v™ 1 ). (A.3) 

In IA.3K m a is the mass of particle a and = r a is its relative velocity with respect to the center of 
mass. Since the position of the center of mass x is given by 



(A.4) 



E m a 

we have the following identities which will be used later: 

J2 m c.r a = 0, ^m a < 8l =0. (A.5) 

a a 

Next, we apply a tensor product with r a to both sides of jA.2t to give 

r a ® fa =ra®p a - (A.6) 



On using the identity 
equation )A.6t becomes 



7(i , Q 0p t> ) = »" 1 0p Q +r„0p a , (A.7) 
at 

^r(ra®Pa) = Wa+2Ta, (A.8) 

at 

where W Q = r a ® fa is the virial tensor and T a = 5("q c1 ® Pa) is the kinetic tensor of particle 
a. Equation )A.8l is called the dynamical tensor virial theorem. This "theorem", which is simply an 
alternative form for the balance of linear momentum, becomes useful in a modified form after making 
the assumption that the atoms in the system are in stationary motion. This means that there exists a time 
scale t, which is short relative to macroscopic processes but long relative to the characteristic time of the 
atomic system, over which the atoms remain close to their original positions with bounded positions and 
velocities. This condition is satisfied for a system of atoms undergoing thermal vibrations about their mean 
positions as expected in a solid at moderate temperature. To exploit this property of the system, we define 
the time average of any quantity / over the time r as 



- r /w* 

T JO 



} = - f(t)dt (A.9) 



and apply this averaging to jA.8t : 

1 



-{r a ® Pa) 



Wa+2Ta- (A. 10) 



Assuming that r a ® Pa is bounded, and a separation of time scales between microscopic and continuum 
processes exists, the term on the left-hand side can be made as small as desired by taking r sufficiently 
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large. Therefore the above equation is reduced to W a = — 27~c« ■ Summing over all particles gives the 
tensor virial theorem: 

W=-2T, (A. 11) 

where W is the time-averaged virial tensor of the system, 

W=^T- a <g>/ a , (A. 12) 



and T is the time-averaged kinetic tensor of the system, 



(A. 13) 



The expression for T can be simplified by substituting in )A.U and noting that due to separation of time 
scales x is constant on the atomistic time scale t, so that 



(A. 14) 



2 

The term in the square brackets is zero due to jA.5t . and thus 

T = - ^2m a v™ 1 ® w™ 1 . (A. 15) 



It is important to note that the virial theorem in )A. 1 It applies equally to continuum systems at rest as 
well as those that are not in macroscopic equilibrium and are therefore in a state of motion. This statement 
hinges on the separation of scales assumption according to which continuum motion occurs so slowly rel- 
ative to atomistic processes so as to be essentially constant on that time scale. 

The virial theorem leads to a definition for stress by considering the idea that the forces on the particles 
in the system can be divided into two parts: an internal part, / ™ , resulting from the interaction of particle 
a with other particles in the system and an external part, /„ xt , due to its interaction with atoms outside 
the system and due to external fields, 

fa =/i nt + /S Xt - (A. 16) 

In terms of continuum variables, the external part of the force can be identified with the traction, t, that the 
surrounding medium applies to the system of particles and the body force, pb, acting on it, where p is the 
mass density and b is body force per unit mass0 Substituting iA. 1 6t into )A.12t . divides the virial tensor 
for the system into internal and external parts, 

W = Wi„t + Woxt = T <* ® + Va ® ^ Xt - (A. 17) 

Rewriting the external virial aQ 

Wcxt = y^r a ® fgx* := / x(g>pbdV+ x®tdA, (A.18) 
„ JQ J8S7 

where f! is the domain occupied by the system of particles and d Q is a continuous closed surface bounding 
the particles and separating them from surrounding particles. The variable a; is a position vector within 



43 These density, traction and body forces are pointwise continuum fields, i.e., they are defined at all 
points but do not include the spatial averaging implicit in the macroscopic continuum description. See 
Section l3~T1 for more details on pointwise fields. 

44 We accept this step as an ansatz due to the many assumptions involved, one of which being that 17 is 
large enough to express the external forces acting on Q in the form of the continuum traction t. 
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Q and on the surface df2. Substituting the Cauchy relation, t = crn, where er is the pointwise Cauchy 
stress, into the above equation and applying the divergence theorem, we have 



Wext = / [x®pb + Aiv x (x ®a)] dV = / | er T + x ® (div x a + pb) \ dV. (A.19) 
J (2 



J <> 



We assume that the pointwise fields satisfy the same balance laws as the macroscopic continuum fields. 
Therefore, the term (div^ er + pb) is zero under equilibrium conditions (see 13. 231 ). We therefore have 
that 

Wext = V(T T , (A.20) 
where we have defined the continuum stress field as the average value of u over the domain S2: 



V la 



crdV, 



(A.21) 



Here V is the volume of Substituting )A.20t into iA. 171 and then into the virial theorem in )A. lit 
gives 

<Tav = " ^ [W int + 2T] T . (A.22) 



Substituting in 4A. 15t and the definition of Wi„t from 1A.17K we have 



(A.23) 



This is the virial stress tensor. It is an expression for the Cauchy stress tensor given entirely in terms of 
atomistic quantities. Finally to demonstrate the symmetry of the virial stress, we rewrite /^ nt as the sum 
over its decomposition: 

ft* = D f°p> (A - 24) 

a 

where f a a are the terms in the central-force decomposition corresponding to a potential energy extension 
as explained in Section l3~4l Substituting jA.24t into jA.23K we obtain 



(A.25) 



Now recalling that f a p = — fp a , we have the following identity 

J"! fa/3 <& r a = \^2 (f a $ ® Va + fP a ® V p) = 9 Y f a P ® ( - Va ~ 



(A.26) 



a,0 



Substituting IA.261 into IA.25K we have 



(A.27) 



This expression shows that the virial stress is symmetric, since f a is parallel to r a — ra. 



45 The definition of this volume is somewhat arbitrary. One can possibly define V as the total volume 
of the Voronoi cells of the atoms lying within Q. 
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B Distance geometry 



In Section \3A\ we saw that the interatomic potential energy, Vi nt , of a set of N particles can be defined 
as a function on a 3./V — 6 dimensional manifold embedded in IR^t^^ 1 )/ 2 , where each point on this 
manifold represents an N(N — l)/2-tuple of real numbers which correspond to the distances between 
the iV particles in R 3 . We also noted that these N(N — l)/2 distances must satisfy certain geometric 
constraints in order to be physically meaningful (see footnote [T9V In this appendix, we discuss the nature 
of these geometric constraints. 

An entire field of geometry, referred to as distance geometry, has emerged to describe the geometry 
of sets of points in terms of the distances between them. The subject was first treated systematically by 
Karl Menger in the early twentieth century and summarized in the book by Blumenthal (5). More recent 
references include |4 10 11 25 56]. Our discussion below is partly based on |50|, which provided a 
particularly clear explanation of the subject. 

The key function in distance geometry is the Cayley-Menger determinant, \ : R-^-^" 1 )/ 2 — ► R, 
which is defined as 



X(Cl2, ■ • -,ClJV,C23, C(JV-I)JV) = dct 
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where 



■ r 2 



(B.l) 



(B.2) 



For a system of N points, {xi, . . . , xpf}, embedded in R 3 , the Cayley-Menger determinant evalu- 
ated at the point, (ri2) ■ • • , '"(JV— l)Jv))> where r a p = \\x a — is related to the volume Vn—1 of a 
simplex of N points in an (JV — l)-dimensional space through the relation: 



X(ri2, 



■ !^(JV-l)iV 



(-i) N 



V N _ 1 (x 1 ,. . . , x N ). 



(B.3) 



For TV : 



X (r 12 ) = 2L 2 , 

where L = v/si2 is the length of the segment defined by the two points. For N = 3, 

X(ri2, ri3, r23) = -16A 2 . 
where A is the area of the triangle defined by the three points. For N = 4, 

X(n 2 ,...,r-34) = 288V 2 , 
where V is the volume of the tetrahedron defined by the four points. For N > 5, we must have 

X(7l2. ■ ■ •,r( JV _i) A r) = 0, 

for points in R 3 since any simplex with five or more points has zero volume in three-dimensional spaceF^I 
We now seek to go in the opposite direction. Rather than computing the squared distances {s a p} 
from a set of points and using the Cayley-Menger determinants to compute volumes, we seek to verify that 
a set of distances actually corresponds to a set of points in three-dimensional space. In technical terms, we 
want the points associated with the distances to be embeddable in R 3 . In order for this to be the case the 
following conditions must be satisfied Fl 



(B.4) 



(B.5) 



(B.6) 



(B.7) 



This is easier to visualize in two-dimensional space. In that case, a simplex with four vertices (a 
tetrahedron) would have zero volume since its vertices would be confined to a plane. The same applies to 
higher-order simplexes and the corresponding higher-order volumes. 

47 Actually, a somewhat stronger theorem can be proved. If the points are numbered in such a way that 
the first four points satisfy conditions [T]and|2] then conditions [3]and|4]need only be applied to groups of 
points that include these four points. See [3j for details. 
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Fig. B.l A cluster of 5 particles. The cluster is shown projected onto the plane normal to the plane defined 
by particles 1 , 2 and 3 shown as a horizontal line. Atom 4 lies above this plane. There are two possible po- 
sitions for atom 5 (at the same height above or below the 1-2-3 plane), where the distances {r±2 , ^35} 
are all the same and only distance r.45 differs. The alternative position of atom 5 and corresponding dis- 
tance r45 are shown in parentheses. Based on Fig. 2 in 1391 . 



1. X(r a p, ) < 0, Vet < P < 7, 

2- x( r a/3, r aJ , r aS ,..., r 7i ) > 0, Va < /3 < 7 < 8, 

3- x(,r a 0,r a y,r a s,r ae , . . . ,r Se ) = 0, Va</3<7<<5<e, 

4. x(r a p,r a - l ,r a 5 > r ae ,r ci <; > ...,r < ,£) = > Va</3<7<<5<e<C- 

For example, Condition^states that the Cayley-Menger determinants computed for all distinct triplets of 
points must be negative or zero. Conditions 12141 apply similarly to sets of four, five and six points. The 
above conditions are easy to understand in terms of )B.4t - fB~7l . They enforce the correct sign of areas and 
volumes in three-dimensional space and the degeneracy condition in IB. 71 . 

As an example, let us see how this is applied for a cluster of N = 5 particles. There are (5 X 4)/2 = 
10 distances, and only 3x5 — 6 = 9 degrees of freedom. There is therefore one more interatomic distance 
than is needed to describe the configuration and indeed there is one Cayley-Menger determinant coming 
from condition [3] 

" S12 S13 S14 SIS 1" 
S12 S23 S24 S25 1 

X\Ti2, ■ ■ ■ , T45) = det = 0. 

V ' S14 S24 S34 S45 1 

S16 ^25 «35 *45 1 
.1 1 1 1 10. 

This expression can be expanded out leading to the following explicit expression 139 iW 



5 5 5 5 5 



a ps lS s St s el - 12(s a p) 2 s jS s Si + 3(s a/3 ) 2 (s 7l5 ) 2 ] = 0, (B.8 



a=l /3=1 7=1 5=1 e=l 



where in the above equation s a a = sa a whenever a > /3. For a given set of nine squared distances, say 
{sil , . . . , S35}, )B.8t provides a quadratic equation for the tenth squared distance, S45, 



A(s 45 ) 2 + SS45 + C = 0, 



(B.9) 



8 Note that Martin 1391 has a small typographical error in his relation. 
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where A, B and C are functions of the other squared distances. (See ['30 ] for explicit expressions for A, B 
and C.) This means that there are two possible solutions for S45 (and hence also for r^) when the other 
distances are set. This situation is demonstrated in Fig. IB. II 

The above example shows that although only 9 degrees of freedom are necessary to characterize a 
cluster of 5 particles, selecting a subset of 9 distances is not sufficient. For this reason, any interatomic 
potential energy extension (see Section l3~4l must be expressed as a function of N(N — l)/2 arguments 
and not just an arbitrary subset of 3N — 6 of them. 



C Useful lemmas 

In his 1955 paper on the "Derivation of the fundamental equations of continuum thermodynamics from 
statistical mechanics" 1481 . Walter Noll proves two lemmas that play an important role in his derivation. 
In this section, for completeness, we derive Noll's first lemma and then extend it to arbitrary curved "paths 
of interaction". We then derive Noll's second lemma for this more generalized case. 

Let f(v, w) be a tensor-valued function of two vectors v and w, which satisfies the following three 
conditions: 

1. f(v, w) is defined for all v and w and is continuously differentiable . 

2. There exists a S > 0, such that the auxiliary function g(v, w), defined through 

g{v, w) := f(v, m)\\vf+ s \\wf+ s , (C.l) 

and its gradients V„g and V^g are bounded. 

3. f(v, w) is antisymmetric, i.e., 

f(v,w) = -f(w, v). (C.l) 



Lemma C.l Under the conditions mentioned above, the following equation holdsF^\ 
I f(^,y) dy = — - diva! / / f(x + sz, x — (1 — s)z) ds 

JySm 3 2 7iGR3 Us = 



) zdz. 



(C.3) 



Proof Conditions (1) and (2) guarantee absolute convergence and allow the order of integration to be 
swapped. From )C.2t we have 



/ f(x,y)dy = - f(y,x)dy 

Jy£K 3 Jyem 3 



(C.4) 



Introduce new integration variables: on the left replace y with z = x — y and on the right replace y with 
z = y — x. Thus, 

/ f(x,x-z)dz = - f(x + z,x)dz. (C.5) 

Jz£R 3 JzeR 3 
Note that a minus sign on the left-hand side is dropped by formally reversing the integration bounds. The 
two terms in )C.5t are equal to J y ^ m 3 f(x, y)dy, we can therefore write 



iyeM 

Next, from the chain rule, we have 



f f(x,y)dy = ^-[ [f(x,x-z)-f(x + z,x)]dz. 
Jvem 3 2 Jzcm 3 



(C.6) 



Vi/(a: + sz, x - (1 - s)z) = V„/ + V™/, 



where s£K. Similarly, 



(C.7) 



— f(x + sz,x-(l- s)z) = (V„/ + V™/) z. 

ds 



(C.8 



45 The expression in Noll's paper appears transposed relative to 4C.3t . This is because the gradient and 
divergence operations used by Noll are the transpose of our definitions. See end of Section^ 
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Fig. C.l The property of a path of interaction as mentioned in the definition is illustrated in this figure. 
Two paths of interaction for the pairs (u, v) and (u, v) are shown in the figure. 



Combining flC.7t and IC.81 . we have 



[V«b/(sb + sz, x - (1 - s)z)] z = — f(x + sz,x - (1 - s)z). (C.9) 

(is 



z = /(a: + z,a;) - f(x,x-z). (CIO) 



Integrating the above equation with respect to s from to 1 gives 
Vaj / f(x + sz, a; — (1 — s)z)ds 

Js=0 

Substituting )C. 10) into IC.61 . and using the identity 

(V x T)a = div a! (T(g)a), (C.ll) 
where T(x) is a tensor of any order and a is a vector which is not a function of x, gives iC.3\ . □ 
Lemma lQll provides us a solution F(x), to the equation 

f(x,y)dy = div !c F(x). (C.12) 



/, 



e«3 



It is clear that F(x) is only a single solution out of an infinite set of solutions that differ from it by a 
divergenceless second-order tensor field. The following lemma extends Lemma lC~Tl to accommodate other 
possible solutions of which the solution given in {C3\ is a special case. Before that we give the following 
definition: 

Definition 3 For any u and v £ R 3 , the "the path of interaction" of u and v is a contour that joins u 
and v such that for any u and v in M 3 , where \\u — v\\ = \\u — v\\, the path of interaction of pairs (u, v) 
and (u, v) are related by a rigid body transformation {Q \ c), for some Q 6 SO(3) and c £ R 3 , with 
Q(u — v) = u — V and Q = I, whenever u — v = ii — v. See Fie. \C.l\ 

Basically, this definition enforces the condition that the "path of interaction" is a contour whose shape is 
only a function of the distance between the points that it connects. The shape of the contour for a given 
distance is assumed to be dictated by the nature of bonding in the material. Here we assume this shape to 
be known. (See also footnote|28]on page l351 . 

Let (ei , e.^ , e$ ) be a basis of R 3 . For every £ > 0, let T; : [0, 1] — > R 3 be a continuously differen- 
tiable contour in R 3 such that 

Ti(s)-ei=sl, 0<s<l, (C.13a) 

r ; (o) = (0,0,0); r i (i) = G,o,o). (c.i3b) 

Fig. lC.2l describes the properties of the contour Tj mentioned above. 



By the definition of the path of interaction, it is clear that the shape of any path of interaction can be 
described by the contour Tj . Moreover, from its properties it is possible to explicitly define the path of 



80 



3 3 




i i 

(a) (b) 

Fig. C.2 Frame (a) shows an admissible contour and frame (b) shows an inadmissible contour which 
violates property )C.13al of Yi . 




Fig. C.3 The contour r(.; s, x, z) passes through x and the vector z is the difference between the two 
end points of the contour. It is related to the contour such that any point Yn z n (s) is mapped to 

r(s; s, x, z) (shown above as a hollow dot) byarigid body motion (Qj |c), where represents rotation 
and c = x — Q z Yu z u (s), represents translation. 



interaction as a function of Yj, a given point x through which it passes, and a vector z which connects its 
end points. More precisely, for any < s < 1, x g R 3 and z e R 3 , let r(- ; s, x, z) : [0, 1] — > R 3 
denote a path of interaction in R 3 , such that 

r(s; s,x,z) = x + Q z [r M (s) - T M (s)] , (C.14) 

for some Q 2 g SO (3), satisfying 

0*ei=— n ^ r , (C.15) 

\\z\\ 

as shown in Fig. IC.3l Let us now verify that r qualifies as a path of interaction. It is clear from )C.14t that 
the contours r (s; s, x, z) and (s) are related through a rigid body transformation, such that 

r(s;s,x,z) = x, (C.16a) 

r(l; s, x, z) - T(0; s, x, z) = ~z. (C.16b) 
Fig. lC.3l describes the properties mentioned above. 

From )C.16al . it follows that r passes through the point x and from IC.16M it is clear that z is the 
vector joining its endpoints. Moreover Q z in IC.141 is made independent of s to ensure that the condition 
= J whenever x — y = u — v, in the definition of the path of interaction, is satisfied. 
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For any i£R 3 , let 

C x = {r(.; s, x, z) -. < s < 1,2 £ R 3 }, (C.17) 

denote the set of all paths of interaction that pass through x . Now for every y on a contour r 6 Cm , let y± 
denote the projection of y on the line joining the end points of the path. It is easy to see that y± (s; s, x, z) 
is given by 

y±(s;s, x,z) = x + Q z [s||z||ei - 7|| 2 ||(s)] . (C.18) 
Therefore, from )C.16at we have 



x ± (s, x, z) = y±(s; s,x,z) = x + Q z [s||z||ei - (s) 



(C.19) 



Using )C.15t this simplifies to 

x±(s, x, z) = x — sz — QzY\\ z \\ (s). (C.20) 
We now generalize Noll's Lemma lC.ll to arbitrary paths of interaction. 

Lemma C.2 For given paths of interaction on R 3 X R 3 , and under conditions \l\3\ on f(v, w) given at 
the start of this appendix, the following equation holds: 



Jyen 3 



x, y) dy = i div x f f f(x ± +sz,x ± ~(l~s)z) 

2 J 2£IR 3 Us=0 



where x± (s, x, z) is given by iC.20t . 
Proof From <C~6l we have 



>Q x r( M (s) dsdz, 
(C.21) 



/ f(x,y)dy = \[ [f(x,x- z)- f(x + z,x)]dz. (C.22) 

Jy£R 3 2 ./z£IR 3 



Next, from the chain rule, we have 

Vi/(iei + sz, x± — (1 - s)z) = (V„/ + SJ w f)SJ x X A 
where s£l. From IC.201 we have Vx^x = I. Therefore 

Vaj/(a3x + sz, ajj_ — (1 - s)z) = V„/ + V^/. 

Similarly, 

-J{x ± + sz, x± - (1 - s)z) = (V„/ + V„/) 



From )C.2()t . we have 



dx± 
ds 



dx± 
ds 



-z - QzTLJs). 



Substituting jC.26t into jC.25t . we have 

+ sz, x x - (1 - s)z) = - (V„/ + V w f) Q*T' {s). 

ds 1 1 



(C.23) 



(C.24) 



(C.25) 



(C.26) 



(C.27) 



Combining jC.24t and jC.27t . we have 

-|/(£c ± + Sz, a;x - (1 - s)z) = - [V*/ (a5 X + sz, sj. - (1 - s)z)] Q^^y (5). (C.28) 

Integrating both sides over the interval s £ [0, 1], we have 

f(x±(l, x, z) + z,x±(l,x, z)) — f(x±(0, x, z), x±(0, x, z) — z) = 



r iv./( 

J 3 = 



x± + sz, x± — (1 — s)z)QzTC!,« (s) ds. 



(C.29) 
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Using <C20l . (CA5\ , and property )C.13b> of T t , we have 

x± (0, x, z) = x ; x±(l, x, z) = x. (C.30) 



Substituting, IC.301 into )C.29t and using the identity {CI It . we have 
f(x + z,x) — f(x, x — z) = — div x / f{x±+sz, xj_ — (1 — s)z)®Q z T^ z ^(s)ds. (C.31) 

J s — 

If we now substitute the above equation into )C.22t , the lemma is proved. □ 

Lemma C.3 Let Si C K 3 , with a piecewise smooth surface S. For given paths of interaction on M 3 X R 3 , 
and under the conditions \l\3\ on f(v, w) given at the start of this appendix, the following equation holds: 



/ / f{x,y)dydx 



ten J yes 

Iff f f{x ± +sz,x x -(l-s)z)(Q z rL l ,(s)-n)dsdzdS(x), (C.32) 



where x± (s, x, z) is given by hC.20\ and Sl c := 
Proof We immediately see that due to antisymmetry of f(v, w), we have 

f f f(x,y)dydx = 0. (C.33) 
JxeQ JyeQ 



cgfi Jye 

Therefore, we have 



/ / f(x,y)dydx= f f f(x,y)dydx. (C.34) 

J seen J yen c J xes~i J yen 3 

From Lemma lQ2l we have 

/ f(x,y)dy = diva,g(x), (C.35) 

JyefL 3 



with 



1 



g (x) = - I f f{x± + sz,x x -{l-s)z)®Q z TL ll (s)dsdz, (C.36) 
where x x (s, x, z) is given by dC.20t . From the divergence theorem we have 

/ div x g(x) dx = / g(x) ■ n{x) dS (x) . (C.37) 
Jmen Js 

Using (CM\ - (C3T\ , we obtain fC32i . □ 

D Derivation of the Hardy stress from the Murdoch-Hardy procedure 

The Hardy stress, cr^ v , was derived in Section l4~2l subiect to the condition: 

div^ er^ iV (sc,t) = ^ f al3 w(x a - x). (D.l) 

a, P 

Recall that the expression for er^ v given in )4.32t was obtained by using the generator function g H 
(see 14.291 ). which is actually a distribution. Due to the inherent obstacles present in this procedure (see 
footnote [34] on page 1471 . we now derive the expression for the Hardy stress given in I4.32K by avoiding 
the use of the Dirac delta distribution. 



S3 



Definition D.l Consider the following definitions for the function rj(x) and associated functions taken 
from d): 
1. Define a molhiier r) e C*°°(R 3 ) by 



(r):= (c e x P ( F ^ T )if|H|<i, 

I if ||t-|| > 1, 

where the constant C > is selected so that J r3 rj dr = 1. 
2. For each e > 0, set 

7? e (r) := ± v (^). (D.3) 
The family of functions r) t are C°° and satisfy 



dr = 1. (D.4) 



The support of t) f is contained in a ball of radius e centered at 0. 
3. If the function h : R 3 — > R is locally integrable, define its mollification h t (r) as 

h c (r):=f Ve (r-y)h{y)dy. (D.5) 

We use the following property of mollifiers in our derivation: 

h c —> h almost everywhere as e — > 0. (D.6) 
For the proof of flD.6t . refer to 1181 . Equation jD.lt can be rewritten as 



diVxd-H v^i*) = ^2 fapV w ( x a ~ x)w(x a - x) , (D.7) 

since w(x a — x) > 0. Now, since \Jw(xg — j/) is locally integrable, using property ID.6K we have 



y/w{x ol - x) = lim / Jw(xg — y)rj e (xB — x a + x — y)dy, (D.8) 

where referring to ID. 51 . h(y) = J w(xp — y) and r + x. Using ED, can be 

rewritten as 

=: lim J2 f g e afj (x,y,t)dy. (D.9) 



Now, note that the function is anti-symmetric with respect to its arguments, x and y, for each e > 
and satisfies all the necessary conditions for the application of Lemma lcTl in Appendix|C] Therefore, from 
IC3l it follows that 

diva, v (x, t) = lim (-^diva; / / f a 0\/w(x a - x - sz)w{xa - x + (1 - s)z) 



z dz 



X TJe (xp — X a + z) ds 

= lim V] f f ( f a 0\/w(x a - x - sz)w(x l3 - x + (1 - s)z) 

e_>0 ^ V ^ ji 3 J s =o V v 

X rj e {xp — x a + z)z ds dz I , (D.10) 



84 



where in the last equality we have used the identity given in IC.1U and the fact that ry E is independent of 
x in the above equation. Since w is positive, and if w ,s a continuously differentiable function, it follows 
that 

Vi (^fal3^w(Xa —X — SZ)W(X0 - X + (1 - s)z)J Z (D.ll) 

is a locally integrable function of z. Therefore, by property |P.6l and using the property, t) t (r) = r] e (— r), 
we have 

1 f 1 

div a! d-^ lv (x,t) = -- V, / Vj, [/^^(ajc - x - s(x a — xp))] (x a - xp) ds 

1 

nZl I -[f a Blv((l - s)Xu+ SXp - X)%) (X a - Xg)]d,S 

2 ■'^=0 



diva. 



2 

L a„ 



(D.12) 



where referring to JD.5I . r = x a — xp, y = z, and h(z) is given in JD.1U . Comparing both sides of 
(D. 12K we arrive at the expression for the Hardy stress tensor given in )4.32t . 
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